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Abstract 

The  quantification  of  uncertainties  is  critical  when  systems  are  nonlinear  and  have  uncertain  terms  in  their 
governing  equations  or  are  constrained  by  limited  knowledge  of  initial  and  boundary  conditions.  Such 
situations  are  common  in  multiscale,  intermittent  and  non-homogeneous  fluid  and  ocean  flows,  and  other 
non-linear  dynamical  systems.  The  Dynamically  Orthogonal  (DO)  field  equations  provide  an  efficient  time- 
dependent  adaptive  methodology  to  predict  the  probability  density  functions  of  such  dynamics.  The  present 
work  derives  efficient  computational  schemes  for  the  DO  methodology  applied  to  unsteady  stochastic  Navier- 
Stokes  and  Boussinesq  equations,  and  illustrates  and  studies  the  numerical  aspects  of  these  schemes.  Semi- 
implicit  projection  methods  are  developed  for  the  mean  and  for  the  orthonormal  modes  that  define  a  basis 
for  the  evolving  DO  subspace,  and  time-marching  schemes  of  first  to  fourth  order  are  used  for  the  stochastic 
coefficients.  Conservative  second-order  finite-volumes  are  employed  in  physical  space  with  new  advection 
schemes  based  on  Total  Variation  Diminishing  methods.  Other  results  specific  to  the  DO  equations  include: 
(i)  the  definition  of  pseudo-stochastic  pressures  to  obtain  a  number  of  pressure  equations  that  is  linear  in 
the  subspace  size  instead  of  quadratic;  (ii)  symmetric  advection  schemes  for  the  stochastic  velocities;  (iii) 
the  use  of  generalized  inversion  to  deal  with  singular  subspace  covariances  or  deterministic  modes;  and  (iv) 
schemes  to  maintain  orthonormal  modes  at  the  numerical  level.  While  (i)  and  (ii)  are  specific  to  fluid  flows, 
(iii)  and  (iv)  are  important  for  any  system  of  equations  discretized  using  the  DO  methodology.  To  verify 
the  correctness  of  our  implementation  and  study  the  properties  of  our  schemes  and  their  variations,  a  set 
of  stochastic  flow  benchmarks  are  defined  including  asymmetric  Dirac  and  symmetric  lock-exchange  flows, 
lid-driven  cavity  flows,  and  flows  past  objects  in  a  confined  channel.  Different  Reynolds  number  and  Grashof 
number  regimes  are  employed  to  illustrate  robustness.  Optimal  convergence  under  both  time  and  space 
refinements  is  shown  as  well  as  the  convergence  of  the  probability  density  functions  with  the  number  of 
stochastic  realizations. 

Keywords:  Uncertainty  Quantification,  Navier-Stokes,  Boussinesq,  Dynamical  Orthogonality,  Projection 
Methods,  Total  Variation  Diminishing,  Error  Subspace  Statistical  Estimation,  Ocean  Modeling,  Data 
Assimilation 


1.  Introduction 

Quantifying  uncertainty  is  becoming  increasingly  important  in  many  scientific  and  engi¬ 
neering  applications.  This  is  in  part  because  the  accuracy  of  an  answer  is  now  often  as  critical 
as  the  answer  itself.  Our  present  motivation  is  uncertainty  prediction  for  computational  fluid 
dynamics  (CFD)  applications,  specifically  in  the  context  of  realistic  ocean  predictions.  In 
ocean  dynamics,  it  is  challenging  to  model  multi-scale,  intermittent,  non-stationary  and 
non-homogeneous  uncertainties.  Already  a  single  evaluation  of  an  ocean  model  is  costly  and 
straightforward  stochastic  modeling  methods  are  prohibitively  expensive  [36,  46],  particularly 
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when  dealing  with  longer-term  unsteady  nonlinear  dynamics.  Fortunately,  the  recently  de¬ 
veloped  Dynamically  Orthogonal  (DO)  held  equations  [50,  49,  51]  provide  efficient,  tractable 
equations  for  uncertainty  prediction  in  large-scale  CFD  and  ocean  applications.  While  these 
DO  equations  have  been  solved  numerically  using  a  simple  finite- difference  scheme,  the  spe¬ 
cific  properties  of  the  DO  equations  warrant  novel  integration  and  discretization  schemes. 
Hence,  our  present  goals  are  to  derive  efficient  computational  schemes  for  the  DO  methodol¬ 
ogy  applied  to  unsteady  stochastic  Navier-Stokes  and  Boussinesq  dynamics,  and  to  illustrate 
and  study  the  numerical  aspects  of  these  schemes.  While  we  are  specifically  focused  on 
incompressible  fluid  hows,  the  use  of  generalized  inversion  to  deal  with  singular  subspace 
covariances  or  deterministic  modes,  and  schemes  to  maintain  orthonormal  modes  at  the 
numerical  level  are  relevant  to  any  system  solved  using  the  DO  method. 

Stochastic  modeling  approaches  can  be  categorized  as  either  non-intrusive  or  intrusive. 
Non-intrusive  approaches  have  the  advantage  that  the  deterministic  version  of  a  model  can 
be  used  to  generate  an  ensemble  of  sample  solutions  from  which  the  statistics  can  be  cal¬ 
culated.  The  disadvantage  is  that  a  large  number  of  samples  are  often  required,  leading  to 
large  computational  costs.  Intrusive  approaches  require  the  update  of  existing  codes  or  the 
development  of  a  new  code,  where  the  resulting  system  of  equations  can  be  larger  than  the 
original,  deterministic  system.  However,  intrusive  methods  are  usually  computationally  less 
expensive  than  non-intrusive  methods,  and  the  statistics  are  explicitly  available.  Below  we 
give  a  brief  review  of  existing  methods,  focusing  on  their  computational  aspects.  For  more 
complete  references  and  reviews  we  refer  to  e.g.:  Ghanem  and  Spanos  [16],  Kloeden  and 
Platen  [24],  Doucet  et  al.  [9],  Mathelin  et  al.  [43],  Eldred  et  al.  [10],  Jakeman  and  Roberts 
[22],  Najm  [47],  Xiu  [66,  67],  Le  Maitre  and  Knio  [27]. 

The  non-intrusive  Monte-Carlo  method  provides  access  to  the  full  statistics  of  the  prob¬ 
lem.  Its  computational  cost  does  not  strictly  depend  on  the  size  of  the  system,  but  more 
on  the  number  of  truly  independent  random  variables,  and  convergence  rates  are  often  pro¬ 
portional  to  the  square  root  of  the  number  of  samples.  The  efficiency  can  be  improved  for 
example  by  using  more  elaborate  Monte-Carlo  schemes  [e.g.  9],  including  particle  filters  or 
mixtures  of  weighted  kernels,  e.g.  Gaussian  kernels  [5,  45].  Nonetheless,  a  large  number  of 
function  evaluations  are  needed  due  to  the  slow  (square  root)  convergence,  which  can  limit 
accuracy  in  large-scale  applications. 

The  Polynomial  chaos  expansion  (PCE),  pioneered  by  Ghanem  and  Spanos  [16]  and 
based  on  the  theory  by  Wiener  [64,  2,  65],  has  become  popular  because  it  can  represent  and 
propagate  large  uncertainties  through  complex  models.  Both  non-intrusive  [e.g.  63,  25,  21, 10] 
and  intrusive  [e.g.  7,  26,  61,  42,  47]  versions  have  been  employed,  but  both  can  suffer  from 
the  curse  of  dimensionality.  That  is,  the  PCE  of  a  dynamical  model  scales  as  ,  where 
p  is  the  largest  degree  of  the  polynomials  used,  and  s  is  the  number  of  independent  random 
variables.  For  problems  that  require  large  p  (e.g.  non-Gaussian)  and  large  s  (e.g.  large  ocean 
or  fluid  simulations),  the  number  of  function  evaluations  (non-intrusive),  or  the  size  of  the 
system  of  equations  (intrusive),  can  quickly  become  prohibitive.  For  s  larger  than  p,  the 
storage  of  a  PCE  scales  as  0(sp )  while  its  computational  cost  for  Navier-Stokes  flows  scales 
as  0(s2p )  due  to  the  quadratic  nonlinearity  of  advection  terms.  The  large  cost  of  these 
methods  have  prompted  the  use  of  non-Gaussian  random  variables  [48],  the  development  of 
generalized  PCE  [68]  to  speed  up  convergence  in  the  polynomial  degree  (i.e.  reduce  p),  and 
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the  development  of  adaptive  schemes  that  only  evaluate  the  necessary  terms  in  the  PCE  (Li 
and  Ghanem  [38]). 

PCEs  have  been  successful  in  many  CFD  applications.  In  the  case  of  unsteady  incom¬ 
pressible  fluid  dynamics,  Le  Maitre  et  al.  [28]  used  a  PCE  scheme  to  study  mixing  in  a 
two-dimensional  (2D)  microchannel  and  improved  the  efficiency  of  their  solution  scheme 
by  decoupling  the  velocity-pressure  equations  using  a  projection  method.  Wan  and  Karni- 
adakis  [62]  studied  the  long  term  behavior  of  a  generalized  PCE  first  using  the  ID  advection 
equation  and  then  a  2D  noisy  flow  past  a  circular  cylinder.  The  authors  showed  that  multi¬ 
element  generalized  PCEs  can  significantly  improve  the  accuracy  for  long  time  integration, 
but  caution  that  the  cost  for  large  s  remain  high.  Other  applications  include  fluid-structure 
interactions,  [e.g.  69],  turbulence  [e.g.  41],  and  aerodynamics  [e.g.  54],  Other  examples  are 
also  studied  in  the  above  references. 

Motivated  by  the  multi-scale,  intermittent  and  non-homogeneous  uncertain  ocean  fields, 
the  Error  Subspace  Statistical  Estimation  (ESSE)  method  was  developed  [35,  29,  30].  It 
uses  a  Karhunen-Loeve  (KL)  expansion  [23,  39]  but  with  time  varying  and  adaptive  basis 
functions.  This  generalized  KL  expansion  is  initialized  by  a  multi-scale  scheme  [31]  and 
evolved  using  stochastic,  data-assimilative  and  adaptive,  Monte-Carlo  ensemble  schemes. 
The  computational  cost  of  ESSE  predictions  scales  as  the  size  q  of  the  Monte-Carlo  ensemble 
required,  i.e.  O(q).  This  ensemble  size  q  varies  with  time  and  is  linked  to,  but  a  bit  larger 
than,  the  evolving  size  of  the  error  subspace  itself  s  (which  gives  the  storage  cost  scaling  in 
O(s)). 

The  DO  equations  for  dynamically  evolving  stochastic  fields  [50,  49]  were  derived  to 
approximate  the  Fokker-Planck  equation  (or  Liouvillc  equation  if  no  stochastic  forcing  is 
used)  and  capture  the  dominant  stochastic  subspace  while  being  computationally  tractable. 
The  DO  methodology  also  starts  from  a  truncated  generalized  Karhunen-Loeve  expansion 
but  derives  the  governing  equations  for  the  mean,  the  modes  and  their  coefficients.  In  this 
derivation,  a  key  condition  is  imposed:  the  rate-of-change  of  the  stochastic  subspace  is  dy¬ 
namically  orthogonal  to  the  subspace  itself.  The  DO  subspace  basis,  i.e.  the  DO  modes,  as 
well  as  the  probability  density  functions  (pdfs),  i.e.  the  stochastic  coefficients,  thus  evolve 
only  according  to  the  dynamics  of  the  system.  This  renders  the  DO  decomposition  efficient 
and  limits  divergence  issues.  Consequently,  the  computational  scaling  is  only  dependent  on 
the  number  of  random  variables,  s,  even  when  dealing  with  non-Gaussian  processes:  specifi¬ 
cally,  the  storage  scales  as  O(s)  and  computational  cost  as  0(s2)  for  Navier-Stokes  equations. 
The  size  s  is  in  general  a  function  of  time  so  as  to  adapt  to  the  dynamically  evolving  un¬ 
certainties  and  boundary  conditions  [51].  The  DO  methodology  has  been  applied  to  several 
Navier-Stokes  flows  and  their  stochastic  dynamics  has  been  studied,  including  mean-mode 
and  mode-mode  energy  transfers  for  2D  flows  and  heat  transfers  [49,  52],  However,  the  DO 
method  has  not  yet  been  applied  to  Boussinesq  flows,  and  the  numerical  challenges  of  the 
DO  decomposition  for  the  stochastic  Navier-Stokes  and  Boussinesq  equations  have  not  been 
thoroughly  examined  nor  resolved.  This  explains  the  need  for  the  present  study. 

In  what  follows,  the  equations  for  incompressible  stochastic  Navier-Stokes  and  Boussinesq 
dynamics  are  given  (§2)  and  their  DO  decomposition  is  outlined  in  Appendix  A.  In  §3,  the 
discretization  in  time  is  developed,  discussing  explicit  and  implicit  schemes.  Importantly,  we 
first  combine  DO  pressure  terms  and  define  new  “pseudo-stochastic  pressures”  in  the  original 
continuous  DO  equations.  With  this  new  definition,  given  that  the  pressure  Poisson  equa- 
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tions  and  the  corresponding  matrix  inversions  often  dominate  the  computational  cost,  the 
cost  of  DO  integration  schemes  is  substantially  reduced:  instead  of  scaling  as  0(s2)  [51],  it 
scales  as  O(s)  (as  long  as  the  solution  of  the  pressure  dominates  the  cost  of  the  scheme).  The 
time  integration  schemes  are  then  derived.  For  the  mean  and  the  modes,  we  employ  projec¬ 
tion  methods  [17],  outlining  schemes  of  first  and  second  order.  For  the  stochastic  coefficients, 
we  obtain  several  time-marching  schemes  of  first  to  fourth  order,  and  briefly  discuss  their 
extension  to  the  cases  of  additive  and  multiplicative  random  forcing.  The  discretizations  of 
the  physical  space  and  stochastic  subspace  are  given  in  §4.  For  the  former,  the  discretiza¬ 
tion  of  diffusion  operators  is  straightforward.  That  of  advection  operators  requires  special 
attention:  since  the  deterministic  DO  modes  have  arbitrary  signs,  how  to  apply  upwinding 
based  on  total  variation  diminishing  properties  is  a  key  question  we  investigate.  For  the 
stochastic  subspace,  a  number  of  possible  discretizations  are  outlined,  including  the  direct 
Monte-Carlo  scheme.  In  §5,  the  questions  of  how  to  deal  with  singular  covariances  and  how 
to  maintain  orthonormal  modes  in  the  presence  of  numerical  round-off  and  truncation  errors 
are  discussed.  For  the  applications  in  §6,  a  set  of  benchmarks  are  defined  and  utilized  to 
illustrate  the  properties  of  our  DO  numerics  scheme.  Specifically,  a  verification  benchmark 
based  on  an  asymmetric  Dirac-stochastic  lock-exchange  flow  is  defined  and  used  to  test  the 
implementation.  A  symmetric  stochastic  lock-exchange  is  then  employed  to  evaluate  the  new 
advection  schemes  for  DO  modes.  The  spatial  and  temporal  convergence  is  studied  with  a 
stochastic  lid-driven  cavity  flow.  The  discretization  of  the  stochastic  coefficients  is  examined 
using  a  flow  over  a  square  cylinder  in  a  confined  channel.  Each  flow  benchmark  is  purposely 
chosen  to  be  different  in  part  to  illustrate  the  robustness  of  our  DO  numerics.  We  also  expect 
that  such  benchmarks  can  serve  as  standard  tests  for  future  schemes  and  implementations. 
Lastly,  conclusions  and  discussions  are  in  §7. 

We  note  that  a  condensed  version  of  this  report  was  submitted  to  the  Journal  of  Com¬ 
putational  Physics,  and  is  currently  under  peer  review. 

2.  Stochastic  Dynamically  Orthogonal  Boussinesq  Equations 

This  section  defines  the  differential  equations  that  we  solve  and  then  briefly  outlines  the 
stochastic  DO  methodology.  The  temporal  and  spatial  discretizations  are  derived  in  the 
subsequent  section. 

The  deterministic  components  of  the  partial  differential  equations  (PDEs)  that  we  solve 
on  a  domain  V  are  non-dimensional  Boussinesq  equations1,  in  the  same  form  as  in  Hartcl 


lrThe  dimensional  variables,  denoted  with  a  hat,  have  been  non-dimensionalized  using:  t  =  t.y  j^;  x  =  xft.; 

U  =  U  \  / gr  h :  f)  =  /)min  “t“  P  (pmax  Pulin') 
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et  al.  [19], 


V  •  u  =  0,  xeD, 

^-^V2u  =  -V.(uu)-Vp  +  pe»,  xe®,  (1) 

%-^mv2p=-v  (up)’  xe®- 

The  non-dimensional  variables  are:  u(x,  t)  =  [u,  v,  w] ,  the  velocity  in  3D;  p(x,  t),  the  density; 
and,  p(x,  t),  the  pressure.  The  vector  e9  is  a  unit- vector  in  the  direction  of  gravity,  (x,  t)  are 
the  non-dimensional  space  and  time  variables,  Gr  =  is  the  Grashof  number  which  is  the 
ratio  of  buoyancy  forces  to  viscous  forces,  Sc  =  v/K  is  the  Schmidt  number  which  is  the 
ratio  of  kinematic  viscosity  v  to  molecular  diffusivity  K  for  the  density  field,  g'  = 

is  the  reduced  gravity,  and  h  is  the  vertical  length-scale.  In  what  follows,  we  denote  the  total 
dynamical  rate-of-change  in  the  prognostic  eqns.  (1)  for  velocity  and  density  by  Cu  and  £p, 
respectively,  i.e.  ^  =  £u  and  ||  —  £,p. 

The  latter  prognostic  equation  for  density  originates  from  the  thermodynamic  energy 
equation  and  an  equation  of  state  (it  arises  from  another  form  of  the  Boussinesq  approxima¬ 
tion  frequently  used  in  ocean  modeling  which  retains  the  temperature  and  salinity  fields  as 
state  variables,  e.g.  [6,  18]).  We  emphasize  that  for  problems  without  density-driven  flows, 
V Gr  =  Re ,  that  is,  the  square  root  of  the  Grashof  number  is  the  Reynolds  number.  The 
approach  and  numerical  schemes  that  we  derive  in  this  manuscript  are  directly  applicable  to 
the  Navier-Stokes  equations. 

We  are  interested  in  solving  eqns.  (1)  in  their  stochastic  form.  We  thus  introduce  the  set 
of  random  events  cu  belonging  to  a  measurable  sample  space  D  and  consider  the  stochastic 
velocity,  density  and  pressure  fields:  u(x,  t;  cu);  p(x,  t;  cu)  and  p(x,  t;  cu).  This  leads  to  stochas¬ 
tic  dynamical  rates-of-change  £u  and  Cp.  If  these  rate-of-changes  are  themselves  uncertain, 
for  example  due  to  parameter  or  model  uncertainties,  then  they  also  depend  explicitly  on  cu. 
In  this  study,  we  mostly  focus  on  uncertainties  arising  due  to  uncertain  initial  conditions. 
We  define  general  stochastic  initial  conditions  as 


u  (x,  0;  cu)  =  u0  (x;  cu)  ,  xeP,  cu  G  D, 
p  (x,  0;  cu)  =  p0  (x;  cu)  ,  xeP,  cu  G  D, 

and  stochastic  boundary  conditions  as 


u  = 

<9u 

dn 


gD  (x,t;cu) , 
gTv  (x,  t;  cu) , 


P  =  9dp  (x,f;cu) , 

dP  (  +  \ 

w~  —  gNp  (x,  t;  uj)  , 


x  G  dVD, 

co  G 

x  G  0VN, 

co  G 

x  G  dVDp , 

co  G 

x  G  dVNp, 

co  G 

(2) 


(3) 


where  the  boundary  conditions  are  separated  into  Dirichlet  and  Neumann  conditions  for  the 
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velocity  and  density  fields  (pressure  boundary  conditions  are  considered  later).  The  resulting 
multivariate  stochastic  eqns.  (l)-(3)  define  the  problem  to  be  solved.  As  in  the  deterministic 
case,  specifics  of  the  solution  depend  on  the  initial  and  boundary  conditions  chosen. 

The  DO  decomposition  of  these  equations  can  be  obtained  from  [49,  50]  and  a  summary 
is  provided  in  Appendix  A.  In  short,  the  DO  methodology  begins  with  a  generalized 
Karhunen-Loeve  expansion  truncated  to  s(t)  terms  [51].  The  vector  of  prognostic  state 
variables  $(x,  t;  u)  =  [u,  p]T  is  decomposed  into  the  sum  of  a  deterministic  mean  component 
<l>(x,  t),  with  s  deterministic  modes  <&j(x,  t),  each  mode  multiplied  by  a  stochastic  coefficient 
Yi(t',uj).  This  decomposition  is  first  substituted  into  eqns.  (l)-(3).  The  DO  condition,  the 
rate-of-change  of  the  stochastic  subspace  is  dynamically  orthogonal  to  the  subspace  itself,  is 
then  utilized.  Orthogonality  is  defined  by  the  spatial  inner-product  (a,  h)v  =  jv  JA(aW)  (I'D 
for  arbitrary  vectors  of  spatial  functions  a  =  [a1,  a2, . .  .]T  and  b  =  [b1,  b2, . .  ,]T.  In  general, 
we  note  that  this  definition  of  the  inner  product  assumes  that  the  different  components  of 
the  state  vector  have  been  properly  normalized  [33,  52],  This  is  not  guaranteed  from  the 
simple  deterministic  non-dimensionalization  used  in  eqn.  (1).  In  fact,  an  additional  stochastic 
normalization  is  usually  needed,  reflecting  the  stochastic  initial  and  boundary  conditions. 
After  some  manipulation  (see  Appendix  A  and  [50,  49])  time-evolution  equations  for  the 
mean,  modes,  and  stochastic  coefficients,  which  are  completely  determined  by  the  dynamics, 
are  obtained.  A  major  contribution  of  this  manuscript  is  to  derive  efficient  discretizations 
in  time  and  space  for  these  equations  and  to  evaluate  the  resulting  computational  schemes 
through  a  set  of  new  benchmarks  for  stochastic  Boussinesq  dynamics. 

3.  Semi-implicit  Time  Discretization 

Solving  the  deterministic  version  of  the  system  of  equations  (1)  implicitly  in  time  of¬ 
ten  requires  not  only  a  large  matrix  inversion  at  each  time-step,  but  also  iterations  at  each 
time-step  to  deal  with  the  non-linear  advection  terms,  e.g.  [13].  Discretizing  their  stochastic 
version  (l)-(3)  using  a  brute-force  Monte-Carlo  scheme  would  have  similar  costs  per  realiza¬ 
tions,  hence  a  total  cost  equal  to  that  of  the  deterministic  version  but  multiplied  by  the  size 
of  the  ensemble.  If  a  DO  decomposition  is  used,  solving  the  DO  system  (A.5)-(A.12)  implic¬ 
itly  would  require  a  matrix  inversion  (s2  +  s  + 1)  times  larger  than  for  (1)  since  the  mean  and 
the  modes  are  coupled  through  the  pressure  and  non-linear  advection  terms  (the  number  of 
pressure  equations  are:  s2  for  s,  s  for  s  and  1  for  p).  While  it  is  possible  to  solve  such 
systems,  our  goal  here  is  to  discretize  (A.5)-(A.12)  such  that  the  equations  decouple,  result¬ 
ing  in  an  efficient  solution  scheme.  This  section  describes  how  this  decoupling  is  achieved. 
First  we  explain  why  we  treat  some  terms  explicitly  and  others  implicitly.  We  then  define 
new  pseudo-stochastic  pressures  that  substantially  reduce  computationl  costs,  develop  Pro¬ 
jection  methods  for  DO  equations  so  as  to  split  the  velocity  and  pressure  terms,  and  present 
time  marching  schemes  for  the  stochastic  coefficients.  The  complete  time  discretizations, 
combining  these  time  marching  schemes  with  the  projection  schemes,  are  summarized  at  the 
end.  The  spatial  discretizations  of  physical  space  and  of  the  stochastic  subspace  are  given 
in  Sect.  54. 


6 


3.1.  Explicitly  and  Implicitly  Treated  Terms 

Much  of  the  decoupling  is  achieved  by  treating  some  terms  explicitly,  resulting  in  a  semi- 
implicit  scheme.  First,  we  choose  to  advance  the  stochastic  coefficients  explicitly,  because 
then  C y- y:i  and  MYjYmYn  can  be  treated  as  constants  when  evolving  the  mean  and  the  modes, 
and  no  iteration  is  required  to  solve  (A. 5),  (A. 8)  and  (A. 11).  Somewhat  similarly,  we  treat 
the  inner  product  terms  (Qi,Q>j)vQ>j  in  (A. 8)  explicitly  to  avoid  iterations.  Next,  we  treat 
the  non-linear  advection  terms  explicitly,  which  is  often  done  in  the  Projection  method 
community  (e.g.  [17]).  This  does  impose  a  stability  constraint  on  the  time-step  size,  a 
Courant  -  Friedrichs  -  Lewy  (CFL)  condition.  Third,  we  treat  the  linear  diffusion  terms 
implicitly  because  they  do  not  couple  the  equations  and  the  resulting  diagonal-dominant 
matrices  can  be  inverted  efficiently.  While  they  could  also  be  treated  explicitly,  this  imposes  a 
much  harsher  stability  constraint  that  could  result  in  very  small  timesteps.  Thus,  to  partially 
decouple  the  evolution  equations,  we  advance  the  stochastic  coefficients,  inner  product  terms 
and  non-linear  advection  explicitly.  However,  these  equations  are  still  coupled  through  the 
pressure. 

3.2.  The  Pseudo-Stochastic  Pressures 

In  this  section  we  first  review  the  explicit  treatment  of  pressure  which  results  in  a  scheme 
requiring  s2  +  s  +  1  solutions  of  stochastic  Pressure  Poisson  equations  (PPEs)  per  timestep. 
Then,  we  discuss  our  new  definition  of  pseudo-stochastic  pressures  that  reduces  the  expense 
to  s  +  1  and  show  that  it  is  a  valid  definition  that  does  not  change  the  order  of  accuracy. 

One  approach  for  handling  the  pressure,  which  was  adopted  in  Sapsis  and  Lermusiaux 
[50],  is  to  treat  it  explicitly.  This  approach  takes  advantage  of  the  fact  that  the  full  stochastic 
pressure  can  be  recovered  at  any  time  instant  by  taking  the  divergence  of  (A. 5)  and  (A. 8), 
inserting  the  decomposition  (A. 2),  and  using  the  divergence-free  form  on  continuity,  noting 
that  V  ■  =  0.  The  result  gives  the  stochastic  Pressure  Poisson  equations  (PPEs)  for  our 

system 


V2p  =  -V  ■  [V  ■  (uu)  -  pe9}, 

V2pi  =  -V  ■  [V  •  (uu*)  +  V  ■  (ujti)  -  p^e9} ,  (4) 

V2Pij  =  -V  ■  [V  ■  (Uj-Uj)  +  V  •  (u*Uj)]. 

A  disadvantage  of  this  explicit  pressure  approach  is  that  it  is  expensive;  to  recover  the 
full  stochastic  pressure,  1  +  s  +  s2  Poisson  equations  need  to  be  inverted,  and  this  would 
often  be  the  dominating  cost  of  the  scheme.  Oceanic  applications  are  expected  to  require 
s  ~  0(  102  —  103)  [33],  which  would  be  very  expensive.  Another  disadvantage  is  that  the 
velocity  computed  with  an  explicit  scheme  will  not  be  divergence-free  after  each  timestep. 

We  can  reduce  the  number  of  PPEs  to  s  + 1  by  defining  new  pseudo-stochastic  pressures. 
The  purpose  of  the  pressure  in  divergence- free  flows  is  to  enforce  continuity.  In  our  stochas¬ 
tic  equations,  continuity  needs  to  be  satisfied  by  the  mean  and  each  modal  velocity  held 
independently  (we  assume  that  the  divergence- free  continuity  equation  is  exact,  without  any 
errors  in  its  form).  Also,  each  of  these  velocity  fields  only  needs  a  single  scalar  held  in  order 
to  satisfy  the  continuity  constraint.  By  inspection  of  equations  (A. 5)  and  (A. 8),  we  therefore 
define  new  pseudo-stochastic  pressures,  which  are  a  combination  of  the  mean,  linear-,  and 


7 


quadratic- modal  pressures: 


P 
Pi 

With  this  definition,  the  quadratic  modal  pressures  are  eliminated  from  (A. 5)  and  (A. 8). 
Thus,  to  evolve  the  mean  and  modes,  we  no  longer  need  to  solve  for  the  quadratic  pressures. 
However,  substituting  (5)  into  the  equation  for  the  evolution  of  the  stochastic  coefficients 
(A. 11),  we  find  that  the  second  term  on  the  right- hand-side  of  (A. 11), 

Pvnn  T  V  ‘  (unUm),  U;)^  (Y.mYn  Cymyn) 

retains  the  projection  of  the  quadratic  stochastic  pressure  terms  in  the  subspace.  At  first, 
this  would  indicate  that  the  quadratic  modal  pressures  are  still  needed,  but  for  commonly 
used  boundary  conditions,  the  projection  cancels,  i.e.  the  inner  product  (Vp^,  Uj)^  is  zero 
(see  [52]).  The  quadratic  stochastic  pressure  term  in  (A.  11)  can  be  dropped  without  any 
penalty.  Thus,  by  defining  new  pseudo-stochastic  pressures  (5),  we  have  shown  that  we 
reduced  the  number  of  PPEs  from  s2  +  s  +  1  to  the  expected  s  +  1. 

3.3.  Projection  Methods  for  the  Mean  and  Modes 

To  obtain  a  numerically  divergence-free  velocity,  we  use  a  Projection  method.  A  large 
number  of  different  Projection  methods  exist;  for  a  recent  review,  see  [17].  Projection  meth¬ 
ods  are  known  for  excellent  efficiency,  but  the  proper  specification  of  boundary  conditions 
remains  a  long-standing  issue.  While  advances  in  this  area  are  still  being  made  [15,  53],  we 
have  chosen  to  use  the  “incremental  pressure-correction  scheme  in  rotational  form”  proposed 
by  Timmermans  et  al.  [58],  which  has  a  proven  temporal  accuracy  [17].  We  first  summarize 
the  classic  versions  of  the  scheme,  then  adapt  them  for  the  mean  and  modes. 

Classic  Projection  Scheme.  Intermediate  velocities  are  first  solved  for  using  a  first  or 
second  order  time-accurate  discretization.  In  both  cases,  the  contributions  from  known 
velocities  (at  previous  or  intermediate  timesteps)  are  written  as  F(u k*,pk*),  where  the  time 
instant  tk *  determines  the  order  of  the  scheme  in  time.  Specifically,  for  first  order  in  time, 
k*  =  k  —  1,  and 

— V/-1  +  F(ufc“1,p*“1),  xGP, 
go,  x  e  dVD , 
g n,  x  £  dVN. 
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1  -—.o  ~  k 

- - —\72uL 

A  t  jGf 

<9ufc 

dn 


=  P  +  CYiYj  Pij , 

=  Pi  +  Cy^Yj  Pmn- 


This  is  followed  by  the  computation  of  the  pressure  correction,  6,  so  as  to  satisfy  continuity, 

V  •  uk 

At  ’ 

0,  x  G  dVD , 

9p,  x  G  dVN, 

where  gp  is  the  prescribed  pressure  difference  at  the  boundary.  With  this  correction,  the  full 
pressure  and  velocity  fields  can  be  recovered  at  the  next  timestep,  i.e. 

uk  =  -  vek, 

Pk  =  pk~'  +  ek  -  uv  •  uk. 

The  first-order  time-accurate  discretization,  k*  =  k—  1,  uses  the  contributions  from  known 
velocities  (at  old  timesteps)  as  F(ufe_1,  pfc-1).  The  second-order  time  integration  scheme  is 
obtained  by  appropriately  choosing  the  guess  values  (•)fc*  with  tk_x  <  tk*  <  tk,  and  using  a 
higher  order  backwards  differencing  method.  For  additional  properties  of  such  schemes,  we 
refer  to  Timmermans  et  al.  [58],  Guernrond  et  al.  [17]. 

Projection  Scheme  for  the  Mean.  We  evolve  the  mean  fields  modifying  the  classic  pro¬ 
jection  method  to  account  for  the  moments  of  the  stochastic  coefficients  and  of  the  chosen 
explicit  and  implicit  terms  (Sect.  §3.1).  Starting  from  the  PDEs  (A. 5)  for  the  mean,  we 
obtain: 


V26  = 

86  _ 
8  n 
6  = 


-fc 

3 _ V2fifc  = 

uk~l 

At  ~ 

VOr 

At 

-  Cky. 

1  l 

V26k  = 

—V 

At 

uk  = 

ufc  —  i 

{V  •  (uu)}‘’  -  Vpfe’  +  pk*e 


%k*  , 


-  c**  {V  •  K-u,)}fc  , 


f  =  f-1  +  6k  -  i/V  •  u*, 


pk 


At  ScVCPr 


fft-i 

72-k  _  P _ 

At 


Vzpk  = 


-  {v  •  (fl/i)}*' -  C? 


with  deterministic  boundary  conditions: 


86 

x  G  8Vd  , 

uk 

=  Sd,  x-  =  0, 
8n 

8uk 

x  G  8T> N, 

8n 

=  g  n,6  =  gP, 

Pk  =  g  DP, 

x  G  8VDp 

8pk 

8n  ~ 

x  G  8VNp 

(6a) 

(6b) 

(6c) 

(6d) 

(6e) 


(7) 
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In  the  above,  the  time  instant  tk* ,  either  previous  or  intermediate,  determines  the  order 
of  the  scheme  in  time.  Specifically,  for  first  order  in  time,  k*  —  k  —  1;  for  second  order,  we 
refer  to  [58,  17].  We  note  that  a  difference  between  a  classic  projection  scheme  and  the  above 
DO  mean  scheme  is  the  presence  of  the  covariances  and  third  moments  of  the  coefficients 
Yf  s  (see  Sect.  §3.4).  A  related  one  is  the  coupling  between  the  differential  equations  for  the 
mean,  mode  and  coefficients  (see  Sect.  §3.5). 

Projection  Scheme  for  the  Modes.  As  for  the  mean,  the  modes  are  evolved  by  modifying 
the  classic  projection  method  for  the  eqns.  (A. 8).  We  obtain: 


u" 


A  t  y/Gf 


1  ~  i.  u; 

V  uf  = 


fc-i 


At 


{V  •  (Ulu)}fc*  -  {V  •  (uu,)}fc*  -  Vjf  +  p*V 


I  j  *  m  *  n 

-  <Q„  iif* 


C YiYj  My.ymyn  {V  •  (unUm)} 


V2ek  =— V  •  U* 

1  At 

u?  =  u?  -  A tvel 
jfi=jfi-'  +  eki-vV-  uf, 


p\ 


At  ScVGr 


A-1 

72  _  ri 

At 


vx  = 


-  {V  •  (u,p)}fe*  -  {V  •  (uA)}* 


with  boundary  conditions: 


^-YjYmYn  ■  (unPm)} 

-  pf  • 


39- 

uf  =  gi,u(=  0),-^  =  °,  X  G  9Dd, 


<9u( 


dn 


y  —  gi,jv(—  0),  6i  —  9i,p{—  0),  x  G  <9Xhv, 

Pife  =  g*,r>p(=  0),  x  G  dT>Dp, 
Qn-k 

-gj-  =  &,«„(=  0).  IE®*,. 


(8a) 

(8b) 

(8c) 

(8d) 

(8e) 


(9) 


where,  again,  the  scheme  is  hrst  order  for  k*  =  k  —  1.  Since  we  focus  on  the  numerics  of  the 
DO  differential  equations,  we  assume  that  the  stochastic  boundary  forcings  are  null,  i.e.  gf  s 
are  null  in  eqns.  (9).  In  the  present  study,  we  do  not  exemplify  problems  with  stochastic 
forcing,  but  we  touch  on  the  numerical  modifications  required  in  such  cases  (see  §3.4  and 
[32]  for  an  example  of  stochastic  ocean  models).  The  time  evolution  of  the  Yfs  is  discussed 
next. 


3Jh  Time  Integration  Scheme  for  the  Stochastic  Coefficients 

To  integrate  the  Ordinary  Differential  Equations  (ODEs)  (A.  11)  for  the  stochastic  coeffi¬ 
cients,  we  assume  that  all  variables  are  available  at  time  tk- 1  and  that  we  integrate  forward 
to  time  tk- 
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First,  we  consider  the  case  where  the  original  governing  differential  equations  only  contain 
uncertain  initial  conditions,  and  no  stochastic  forcing.  This  case  corresponds  to  the  examples 
and  benchmarks  we  consider  later.  It  allows  the  use  of  classic  time-marching  integration 
algorithms  or,  in  other  words,  appropriate  ODE  solvers.  We  consider  a  given  realization  of 
a  coefficient  Yx  at  time  tk~\-  To  integrate  to  time  tk,  we  first  define  an  approximation  to  the 
time-rate  of  change  of  this  coefficient  at  a  given  instant  tk- 1  <  t  <  tk  as  follows: 


<m 

dt 


Y(t) 


^^=v2um  -  V  •  (umu)  -  V  •  (uum)  -  Wpm  +  Pm.e9,  Ui^  Ym(t) 


+ 


^c^_V2pm  -  V  •  (u mp)  -  V  •  (upm)>  P^J  Ym(t) 


-  (V  •  (unum),  u i)v  Ym(t )  Yn(t)  -  CYmYr 


-  (V  •  (u npum),pi)v  Ym(t)  Yn(t)  -  CYmYn  , 


(10) 


where  (,)fc  indicates  that  the  modal  quantities  are  estimates  of  their  values  at  time  t.  The 
numerically  exact  option  is  to  choose  tk-  =  t,  while  the  cheapest  is  to  take  tk-  =  tk- 1 
since  at  that  time,  all  modal  quantities  are  available  from  the  previous  time-step.  Using 
this  derivative  definition,  we  have  employed  several  time-marching  schemes  of  varying  order 
to  advance  the  U  in  time,  including:  a  low-storage  4th-order-accurate  explicit  Runge-Kutta 
integrator  of  the  form 


y(0)  =  Y{t),  y{0)  =  0, 

PY 

y(m)  =  amy (m"1}  +  At  — 

dt 

y-(m)  _  y(m-l)  +  bmy(m) 

Y(t  + At)  =  r(5), 


Yi™-1)  ,t+Cm.At 


for  me  [1 ...  5], 
for  me  [1 ...  5], 


where  the  coefficients  am,6m,cm  are  given  in  Carpenter  and  Kennedy  [3];  a  2nd-accurate 
explicit  Runge-Kutta  scheme  (Heun’s  version) 


K(°)  =  Y(t)  +  At  d^~ 


Y(t  +  At)  =  Y(t)  +  ^  |  ^ 


Y(t) 

Y 
dt 


Y(t) 


dY 

dt 


y(o) 


and  the  first-order-accurate  explicit  Euler  scheme 


dY 

Y(t  +  At)  =  Y(t)  +  At  — 

dt 


Yt 


For  the  Euler  scheme,  k  —  k  —  1.  For  each  stage  of  the  two  RK  schemes,  is  evaluated 
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using  (10).  If  tk-  =  t,  the  modal  quantities  (the  rates  for  the  Yj’s)  are  advanced  to  inter¬ 
mediate  times  using  the  modes  and  mean  PDEs  above  (which  is  expensive).  If  k~  =  k  —  1, 
the  mean  and  modal  quantities  are  not  advanced,  but  the  Ym(ty s  are  still  updated  at  inter¬ 
mediate  times.  Of  course,  the  formal  order  of  accuracy  of  that  scheme  is  limited  by  these 
mean/modal  terms  kept  constant  as  will  be  shown  later  (§6).  While  CymY  could  be  recalcu¬ 
lated  at  each  time  level,  our  simulations  showed  better  results  if  they  were  kept  at  the  same 
time  k~  as  modal  quantities  (in  that  case,  coefficients  and  subspace  remain  consistent). 

Second,  we  consider  the  case  where  the  original  governing  equations  contain  stochastic 
forcing  in  the  form  of  zero-mean  Wiener  processes,  added  linearly  to  the  deterministic  form 
of  the  PDEs  (1).  These  stochastic  forcings  then  appear  in  the  governing  equation  (A.  11)  for 
the  stochastic  coefficients.  Without  entering  details,  the  above  marching  schemes  are  then 
augmented  with  stochastic  integrals  over  tk  —  tk-i  for  the  contribution  of  these  forcings  [24, 
20].  If  the  stochastic  forcings  are  of  multiplicative  form  in  the  original  governing  equations, 
then  expectations  and  coupled  integrals  of  stochastic  terms  are  in  general  also  in  the  mean 
and  mode  equations  of  Sect.  §3.3.  This  is  not  considered  in  the  present  study. 

3.5.  Complete  Time  Integration  Scheme 

We  now  summarize  the  complete  time-discretization  scheme  from  tk- i  to  tk.  Since  we 
have  decoupled  the  equations  (6),  (8),  and  (10),  the  order  in  which  they  are  solved  is  not 
important.  In  fact,  (6),  (8),  and  (10)  could  be  solved  in  parallel.  Presently,  we  employ  the 
following,  serial  approach: 

1.  Calculate/extrapolate  the  statistics  ( C ym yn ,  M y . ym Yn )  to  the  approximated  times  k*,  k~ 
and  store  for  later  use 

2.  Calculate/extrapolate  the  advection  terms  (Vuu,  VuUj  +  Vu,u,  Vujin,-)  to  the  approx¬ 
imated  times  k *,  k~  and  store  for  later  use 

3.  Advance  the  Yds  using  (10),  and  one  of  the  ODE  solvers  in  §3.4 

4.  Advance  the  mean  u  using  (6) 

5.  Advance  the  modes  Uj  using  (8). 

For  the  modes  and  mean,  we  choose  k*  =  k~  =  k  —  1,  resulting  in  a  first  order  accurate 
scheme  for  the  present  applications  (§6).  For  second  order  accuracy,  we  would  use  k *  as 
defined  in  Timmermans  et  al.  [58],  Guermond  et  al.  [17].  For  the  stochastic  coefficients,  a 
higher  order  ODE  solver  may  be  used  in  step  3  which  may  reduce  the  magnitude  of  the 
integration  error.  However,  since  DO  equations  are  coupled,  the  order  of  accuracy  of  that 
step  is  influenced  by  the  choice  of  k*  and  k~  used  for  the  time-integration  of  the  modes  and 
mean.  Hence,  if  k*  =  k~  =  k  —  1,  the  overall  accuracy  of  the  time-integration  is  in  general 
expected  to  be  first  order. 

4.  Spatial  Physical  and  Stochastic  Subspace  Discretizations 

In  this  section,  we  start  by  describing  the  2D  spatial  discretization  of  (6)  and  (8)  on  a 
structured  grid  (§4.1).  We  employ  a  standard  conservative  finite  volume  discretization  of 
second-order.  Special  treatment  is  needed  for  the  advection  by  the  modes  since  they  are 
basis  vectors  in  the  stochastic  subspace  and  thus  do  not  have  a  preferential  direction.  We 
finally  discuss  the  discretization  of  the  s-dimensional  probability  subspace  in  §4.2. 
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4-1 ■  Spatial  Discretization  of  the  Physical  Space 

The  domain  T>  is  discretized  into  a  finite  number  of  non-overlapping  control  volumes. 
Presently,  we  use  rectangular  control  volumes,  which  form  a  structured  Cartesian  grid  that 
has  uniform  spacing  in  the  x  and  y  directions.  Several  choices  exist  for  the  relative  placement 
of  velocity  and  pressure  control  volumes.  Here  we  choose  to  use  a  standard  staggered  C-grid 
[13,  40],  where  the  u-  and  v- velocity  control  volumes  are  displaced  half  a  grid-cell  in  the  x- 
and  //-directions  relative  to  the  pressure  and  density  control  volumes,  respectively. 


4-1.1.  Diffusion  Operator 

The  diffusion  operator  V2  is  discretized  by  using  central  boundary  fluxes.  That  is,  for 
the  field  <f>  at  the  center  (x,  y)  of  a  control  volume  boundary,  we  have  the  following  boundary 
flux  in  the  x-direction  (similar  for  //-direction) 

df  =  <j>{x  +  jf, y)  -  <j>{x  -  ^f,y) 

9X  {x,y) 

which,  on  a  structured  grid,  is  exactly  equivalent  to  the  second-order  accurate  central  finite- 
difference  scheme.  For  the  advection  operator,  the  simple  second-order  central  flux  is  well- 
known  to  be  unstable  (e.g.  [4]),  and  needs  more  careful  treatment,  as  described  next. 


4-1.2.  Advection  Operator 

For  advection  by  a  velocity  component  u,  we  use  a  standard  Total  Variation  Diminishing 
(TVD)  scheme,  with  an  monotonized  central  (MC)  symmetric  flux  limiter  [60].  The  scheme 
can  be  written  for  a  variable  rj  as: 


F(Vi-h)  =  ui-l 


Vi  +  Vi- 1 


2  2 


u{_i 
1  2 


r)i  ~  Vi- 1 


i-i 


At 


u.i—— 

2  Ax 


^(D-i) 


where  the  MC  slope  limiter  T(r)  is  defined  as 

T(r)  =  max  <f  0,  min 


1  +  r 

nun  |  ,  2,2  r 


and  the  variable  r  as 


D_i  = 
1  2 


1 

n  U;  1  + 

2  V  2 


Uj_  l 
1  2 


(^-1+^-2)  +  Uui_i 


Ua  1 
1  2 


(Vi+1  +  Vi 


Ui-dVi  ~Vi-i) 


(11) 


where  iq_  1  is  used  without  interpolation  for  the  density  advection  while  a  second-order 
central  scheme  is  used  for  the  non-linear  u  and  v  advection.  For  more  on  TVD  schemes,  we 
refer  to  the  excellent  text  by  [37]. 

A  possible  issue  with  using  this  scheme  for  DO  equations  arises  from  the  realization 
that  the  absolute  value  of  the  velocity  |u|,  is  a  function  of  the  full  velocity  u  +  YiUi  which, 
depending  on  the  specific  realization,  may  be  either  positive  or  negative;  in  other  words,  |u| 
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is  positive,  but  stochastic.  Fortunately,  we  never  need  the  full  velocity  to  evolve  the  mean 
and  modes  in  (A. 5)  and  (A. 8)  (see  also  (6)  and  (8)).  In  fact,  in  the  case  of  the  mean  velocity 
u,  its  absolute  value  is  deterministic.  Therefore,  the  advection  of  the  mean  velocity  by  the 
mean  velocity  u  •  Vu,  and  the  advection  of  the  velocity  modes  by  the  mean  velocity  u  •  Vu* 
can  use  the  classic  TVD  method  without  modification.  Advection  of  the  mean  by  the  modes 
u*  ■  Vu  and  of  the  modes  by  the  modes  u;  •  VUj,  however,  need  additional  consideration. 
Similar  statements  apply  for  the  advection  of  the  mean  density  and  of  the  density  modes  by 
either  the  mean  velocity  (classic  scheme  is  fine)  or  by  the  modes  (additional  considerations 
are  needed). 

Here  we  propose  three  arguments  for  three  different  advection  schemes  that  can  be  used 
for  these  “advection  by  the  modes”  terms.  First,  if  we  examine  the  equations  from  the 
perspective  of  the  numerical  scheme  only,  a  preferential  advection  direction  will  be  present. 
In  this  case  we  simply  use  the  TVD  scheme  unmodified.  Next,  we  argue  that,  since  the 
stochastic  coefficients  are  zero  mean,  then  the  probability  of  Y)  <  0  is  equal  to  the  probability 
that  Yi  >  0.  This  suggests  that  Uj  should  not  have  a  preferential  direction  of  propagation,  in 
which  case  a  central  differencing  advection  scheme  (CDS)  could  be  used.  Last,  recognizing 
that  the  CDS  scheme  may  cause  oscillations,  we  still  wish  to  limit  the  flux  in  some  way.  A 
direct  approach,  then,  is  to  use  the  TVD  scheme  in  both  directions,  and  average  the  results. 
That  is,  the  present  sign  of  the  modal  velocity  is  used  first  to  calculate  the  advective  terms, 
then  the  negative  of  the  modal  velocity  is  used,  and  the  two  results  are  averaged.  We  call 
this  the  symmetric  TVD  or  TVD*  scheme.  Note  that  the  TVD*  scheme  is  not  a  true  TVD 
scheme  and  is  thus  not  guaranteed  to  be  oscillation-free  for  all  realizations. 

To  better  understand  the  behavior  of  the  TVD*  scheme,  consider  three  scenarios:  when 
calculating  the  flux  using  the  negative  and  positive  modal  velocities  i)  neither  requires  lim¬ 
ited,  ii)  both  require  full  limiting,  iii)  only  one  requires  full  limiting.  For  i),  both  use  a  CDS 
flux  and  the  average  of  the  two  is  a  CDS  flux.  For  ii),  both  use  their  respective  UW  fluxes, 
which,  when  averaged  together,  again  gives  a  CDS  flux.  This  may  appear  undesirable  but 
it  is  sensible  since  the  velocity  should  be  fully  up-winded  from  both  sides  with  equal  prob¬ 
ability,  resulting  in  no  preferential  direction  of  up- winding,  which  gives  CDS.  For  case  iii) 
we  depart  from  a  CDS  scheme  since  here  only  one  side  uses  an  UW  flux  and  the  other  uses 
CDS.  When  averaged  together,  75%  of  the  total  flux  will  come  from  the  up-winded  side,  and 
only  25%  from  the  CDS  side.  Since  any  combination  of  these  three  scenarios  can  happen  in 
practice,  the  TVD*  scheme  is  a  significant  departure  from  the  CDS  scheme. 

The  three  proposed  schemes  are  tested  in  §6.2,  where  we  find  that  the  TVD*  scheme 
performs  best.  Improving  this  TVD*  is  still  possible,  since  minor  oscillations  can  remain. 
A  proper  flux-limited  advection  scheme  may  be  derived  by  looking  at  the  characteristics  of 
the  hyperbolic  parts  of  the  full  system.  However,  since  the  system  has  (s  +  1)  x  d  equations, 
where  d  is  the  dimension  of  the  problem,  this  analysis  is  left  for  future  research. 

4-2.  Discretization  of  Stochastic  Subspace 

The  stochastic  coefficients  exist  in  an  s-dimensional  space,  which  could  become  large, 
from  0(10)  to  O(103)  based  on  our  experience.  In  most  cases,  there  is  also  no  strict  bound 
on  the  value  that  a  stochastic  coefficient  can  take.  Thus,  evenly  dividing  the  s-dimensional 
space  is  not  feasible.  To  discretize  the  uncertainty  subspace,  other  schemes  are  thus  used. 
They  include:  (a)  non-uniform  discretizations  of  the  subspace,  either  using  structured  or 
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unstructured  grids,  possibly  using  schemes  based  on  finite- volumes  or  finite-elements  [16,  41]; 
(b)  solve  a  discretized  version  of  the  PDEs  for  the  probability  densities  of  the  coupled  s 
coefficients,  e.g.  solve  Fokker-Planck  equations  [50];  (c)  parameterize  the  probability  space, 
either  using  Polynomial  Chaos  [16,  68,  42,  14],  in  our  case  extended  to  time-dependent 
polynomials,  or  other  parameterizations  such  as  Gaussian  mixtures  [44,  56]  or  particle  filters 

[9] ,  and  (d),  use  a  Monte-Carlo  approach  [35,  30]. 

We  have  employed  a  few  of  these  schemes.  Here,  we  only  illustrate  a  Monte-Carlo  scheme: 
for  each  realization  the  time-integration  of  Sect.  §3.4  is  used.  In  general,  the  expected 
error  for  the  mean  and  covariance  is  of  O  where  q  is  the  number  of  samples.  For 

efficient  results,  it  is  important  that  samples  are  generated  in  regions  where  the  probability 
is  relatively  high,  based  on  importance  sampling  [9].  At  initial  time  to,  the  distribution  of 
the  Y,j  is  generated  using  the  specified  initial  probability  density  given  by  eqns.  (2).  Here, 
our  focus  is  on  numerical  schemes  and  we  restrict  ourselves  to  simple  distributions  such  as 
Gaussians  or  Dirac  functions. 

Thus,  we  discretize  Y{  by  generating  q  samples,  and  forming  a  q  x  s  dimensional  matrix 
Yr  i.  Each  row  in  Yri  corresponds  to  one  of  the  samples,  and  each  column  corresponds  to  one 
of  the  modes.  During  the  time-integration  step,  each  sample  of  Yri  is  then  advanced  using 

(10) ,  which  is  done  efficiently.  In  all  cases  computed  so  far,  q  can  be  large,  e.g.~  0(  104  — 105), 
but  still  sufficiently  small  such  that  advancing  (10)  for  every  sample  is  not  the  dominating 
cost  of  the  whole  scheme. 

A  drawback  to  this  Monte-Carlo  approach  is  that  rare  events  will  not  be  captured  unless 
a  very  large  number  of  samples  are  used.  Alternative  methods,  such  as  Gaussian  mixtures 
[55,  56]  or  other  approaches  mentioned  above  can  then  be  used  as  an  alternative. 


5.  Implementation  Details 

In  this  section  we  describe  selected  implementation  details  that  arise  due  to  numerical 
issues.  In  particular,  we  discuss  how  to  deal  with  possibly  poorly  conditioned  covariance 
matrices  in  the  stochastic  subspace,  as  well  as  the  orthonormalization  of  the  modes  and 
decorrelation  of  the  stochastic  coefficients. 


5.1.  Dealing  with  a  Singular  Covariance  Matrix 

The  covariance  matrix  may  be  singular  or  poorly  conditioned  if  one  or  more  of  the 
stochastic  coefficients  have  zero  or  very  small  variance  compared  to  other  modes.  This 
situation  for  example  arises  if  a  system  has  deterministic  initial  conditions,  but  becomes 
uncertain  through  forcing,  boundaries,  parameters,  numerical  uncertainties,  or  other  causes. 
The  initial  covariance  matrix  is  then  simply  zero:  its  inverse  is  not  defined.  Special  treatment 
is  thus  needed  for  such  cases  since  the  inverse  of  the  covariance  matrix  is  required  in  (8). 

Fortunately,  this  problem  is  also  common  in  data  assimilation  where  it  is  resolved  using 
generalized  Moore-Penrose  inversions  [1,  35].  In  the  particular  case  of  the  DO  equations,  the 
inverse  of  the  covariance  matrix  is  multiplied  by  the  third  moments  in  (8), 
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which,  for  most  physical  processes,  goes  to  zero  for  the  eigenvalues  of  Cy*y  that  go  to  zero. 
However,  to  ensure  a  numerically  stable  estimate  of  the  inverse  of  the  coavariance  matrix,  we 
employ  a  generalized  Moore-Penrose  inverse,  which  amounts  to  truncate  the  singular  values 
less  than  a  defined  tolerance  or  to  set  them  to  that  tolerance  (with  the  former,  the  inverse 
of  a  zero  covariance,  i.e.  deterministic  initial  conditions,  is  zero).  This  results  in  a  stable 
numerical  simulation,  as  exemplified  for  the  Lock  Exchange  validation  benchmark  in  §6.1. 


5.2.  Orthonormalization 

The  DO  equations  enforce  orthonormal  modes  and  it  is  important  to  maintain  this  prop¬ 
erty  numerically  when  integrating  over  time.  Analytically,  if  modes  are  orthonormal  initially, 
orthonormality  is  maintained  because 


d  (*t ,  *j)P 

dt 


0, 


where  the  DO  condition  ^<l>j ,  =  0,  Vi ,  j  e  [1,2,  ...,s]  was  used.  At  the  discrete 

level,  this  property  is  maintained,  up  to  truncation  and  round-off  errors.  Even  if  the  modes 
are  orthonormal  at  a  given  time-step,  not  all  integration  schemes  over  the  next  time-step 
will  conserve  the  discrete  orthonormality. 

Let’s  first  consider  the  analytical  integration  from  time  t^i  to  t*..  For  the  modes,  we 
have  1  +  A where  the  exact  time  integral  is  denoted  as  A ^ dt . 

For  the  inner  product  at  tj.,  we  then  have: 


<*? .  *5%  =  ( *D  +  Ai  Ai ,  ■  +  At9** 


at 


V 


-h-1 ,  sDL  +  a 


v 


A  t(  $ 


,fc-l  9*  A  ,  A+2  d*i 


+  AO 


dt  /  v  \  dt  ’  dt  ,  v 


Since  ,  &j)v  =  1 ,  l)v  =  5ij,  the  remaining  terms  sum  to  zero, 

A +A tUk~l  d®3  x  '  a,2/—  * 


1  ’  dt 


v 


dt 


.  ,  .  d^i  d&j  . 

+  At  (  —^A  ,  )  =  0 


D 


dt  ’  dt 


v 


(13) 


Considering  now  a pth-order-accurate  discrete  approximation  of  its  error,  £j,  is  of  0(Atp). 
Assuming  that  the  spatial  inner  product  is  computed  exactly,  this  error  in  the  time  in¬ 
tegration  leads  to  an  error,  Sy,  in  the  inner  product  at  time  k,  i.e.  (<£* ,  <&j)^dlscietc  = 
,  <£>J).p  +  Sy.  To  estimate  the  magnitude  of  this  Sy  made  over  one  time  step,  we  assume 
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an  exactly  orthonormal  inner  product  at  tk-i-  By  discrete  integration  to  f*.,  we  then  have: 


Using  (13), 


<*i ,  *i)i  +  %  =  <*?_1 ,  +  A  t(* 


fc-i  ^ 


At  /  +  £ 

\  *  ’  dt  +  J 


At2/^ 


=  At  (^j-1 ,  £i)v  +  At  ,  £j)v 

+  At2  (st ,  +  Af2  (^S  ,  +  At2  <£, ,  y  D , 

~  2At0(l)0(Atp)  +  2At20(Atp)0(Af1)  +  At20(At2p), 
~  C>(Atp+1). 


Therefore,  the  error  in  the  orthonormality  will  always  be  At  smaller  than  that  of  the  numer¬ 
ical  scheme. 

The  orthornomality  can  be  corrected  indirectly  by  enforcing  that  the  solution  and  the 
error  are  numerically  orthonormal,  ($,; ,  +  Sy  =  StJ1  at  the  end  of  a  time  step.  This 

has  to  be  done  with  care:  the  summation  Yr^l  produces  specihc  realizations,  and  changing 
the  basis  without  modifying  the  coefficients  will  change  the  specihc  realizations.  Because  <f>i 
and  Yrj  are  linked,  various  schemes  for  performing  the  orthonormalization  exist.  They  are 
described  in  Appendix  B. 


6.  Numerical  Applications 

In  this  section  we  present  four  benchmarks  used  to  verify  and  study  the  schemes  described 
above.  To  ensure  that  the  implementation  is  solving  the  desired  equations,  we  compare  the 
stochastic  code  to  a  deterministic  code  (§6.1)  for  a  version  of  the  lock-exchange  dynamical 
problem  [19].  Next,  we  examine  three  advection  schemes  proposed  for  DO  equations  in 
§4.1.2,  using  a  symmetric  version  of  the  lock-exchange  problem  (§6.2).  Then,  we  evaluate 
the  spatial  and  temporal  convergence  using  the  lid-driven  cavity  how  (§6.3).  Finally,  we 
study  the  discretization  of  the  stochastic  coefficients  using  the  how  over  a  square  cylinder  in 
a  confined  channel  (§6.4). 

For  evaluating  errors,  we  use  the  L 2  norm.  At  a  single  time  instance,  the  norm  is 
ll*l|2  =  VWMd  for  the  whole  state  vector,  or  ||0||2  =  \J fv  <j)2dT>  for  a  single  component. 

For  the  convergence  studies  (§6.3)  we  also  integrate  over  time,  using  ||<I>||2  =  \J f, jf'  (d*-  dt 

for  the  state  vector,  and  ||Y)||2  =  J /QT  E  [UU]  dt  for  the  stochastic  coefficients. 
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Figure  1:  Lock-exchange  problem:  an  initial  barrier  separating  light  and  heavy  fluid  is  removed,  and  the  flow 
is  allowed  to  evolve.  Uncertainty  in  our  studies  originates  from  not  knowing  the  initial  density  differences 
between  the  fluids.  This  benchmark  is  used  to  verify  the  correctness  of  the  implementation,  not  the  DO 
methodology. 


6.1.  Lock- Exchange  Verification  Benchmark 

The  purpose  of  this  benchmark  is  to  verify  the  numerical  implementation.  For  exam¬ 
ple,  our  use  of  stochastic  pseudo-pressures  does  not  in  theory  introduce  additional  errors. 
However,  one  needs  to  verify  that  their  numerical  implementation  is  accurate,  solving  the 
correct  equations.  Ideally,  problems  with  analytical  solutions  should  be  used  to  verify  a  code, 
however,  constructing  a  valid  analytical  solution  of  (6),  (8)  and  (10)  with  multiple  stochastic 
modes  and  coefficients  is  neither  trivial,  nor  does  it  lend  itself  to  compact  expressions.  In 
this  section  we  address  this  problem  by  defining  a  numerical  benchmark  and  then  using  it 
to  verify  the  present  DO  code. 

6.1.1.  Lock- Exchange  Verification  Benchmark:  Setup 

We  verify  our  stochastic  DO  code  by  comparing  it  to  a  deterministic  NS  code  which  has 
been  thoroughly  verified  [34,  59].  This  deterministic  code  uses  the  same  second-order  Finite 
volume  scheme  and  first  order  backwards  difference  Projection  method  as  the  stochastic 
code.  While  the  DO  code  is  inherently  more  complicated  with  coupled  equations,  the  major 
differences  are  in  the  advection  scheme  for  the  stochastic  modes  (see  §4.1.2  and  §6.2),  and 
in  the  need  to  also  evolve  the  stochastic  coefficients. 

The  deterministic  code  is  used  non-intrusively  with  a  Monte-Carlo  method  to  generate 
an  ensemble  of  independent  realizations.  These  references  are  then  compared  to  realizations 
from  the  DO  code  which  solves  the  coupled  DO  equations.  The  benchmark  (Fig.  [1])  is 
based  on  the  lock-exchange  problem  [19],  where  uncertainty  is  introduced  by  prescribing 
four  possible  initial  density  differences.  In  other  words,  while  the  exact  difference  between 
the  densities  is  unknown,  it  is  known  that  only  four  possibilities  of  equal  probability  exists: 
i.e.  the  pdf  is  initialized  as  four  discrete  Dirac  delta  function.  Essentially,  we  employ  a  single 
run  of  the  DO  code  to  try  to  replicate  four  independent  deterministic  runs.  While  not  a 
practical  use  of  the  DO  method,  it  is  a  challenging  benchmark  for  verifying  its  numerical 
implementation. 

To  capture  all  of  the  uncertainty,  three  DO  modes  (s  =  3)  are  needed  for  the  stochastic 
simulation.  Four  density  difference  were  chosen  because  three  DO  modes  is  the  minimum 
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Mean  Mode  1  Marginal  PDF 


-0.4  0  0.4 


P 

Figure  2:  Initialization  of  the  lock-exchange  problem  (density  field:  mean,  mode  1  and  marginal  pdf,  modes 
2  and  3).  The  initial  velocity  is  zero,  and  stochasticity  is  introduced  through  the  first  (of  three)  orthonormal 
modes  for  the  density.  The  vertical  length  scale  used  for  non-dimensionalization  is  the  half-height  of  the  chan¬ 
nel  (h  =  1).  Initializing  with  four  stochastic  samples  for  the  first  mode,  Yr> i  =  [—0.18,  —  0.06,  0.04,  0.20]T, 
we  have  four  possible  initial  conditions  corresponding  to  A p  =  [1,  0.84,  0.74,  0.62].  The  two  remaining 
modes  are  initialized  as  described  in  the  text,  and  do  not  introduce  additional  stochasticity. 


number  required  to  have  full  energy  interactions  between  the  mean  and  modes  of  the  DO 
simulation  [50,  49].  Also,  to  fully  verify  the  implementation,  we  need  to  ensure  that  the 
initial  pdf  is  non-symmetric  so  that  the  third  moments  are  non- zero  (in  (10)  for  example). 

General  setup:  The  Schmidt  number  is  kept  constant,  Sc  =  1,  and  we  present  results 
for  Gr  =  1.25  x  106  and  Gr  =  4  x  104,  (although  other  Gr  were  studied).  The  four  density 
differences  of  equal  probability  are  A p  =  [1,  0.84,  0.74,  0.62],  with  the  initial  density  profile 
prescribed  by 

A  p 

p(x,y,t  =  0)  =  —  tanh(2  x/lp), 

where  we  take  lp  =  1/64.  Initially  the  velocity  is  zero  everywhere.  Free-slip  boundary 
conditions  are  used  at  the  domain  boundaries  (Fig.  [1]).  The  domain  is  discretized  using 
Ax  =  Ay  =  1/256,  and  At  =  1/512,  which  is  sufficient  resolution  for  these  Grashoff  numbers 

[19]- 

Mean  initialization:  The  mean  density  profile  (Fig.  [2]-a)  uses  the  hyperbolic  tan  function 
specified  above  with  a  mean  density  difference  A p  =  0.8.  The  mean  pressure  and  velocity 
are  zero  everywhere,  initially. 

Mode  initialization:  The  density  profile  for  the  first  mode  is  the  hyperbolic  tan  profile 
above  (Fig.  [2]-b),  but  normalized.  The  two  remaining  modes  are  arbitrary,  since  they  do 
not  introduce  initial  uncertainty  (see  below).  They  are  set  to: 


p2(x,t  =  0) 


p3(x,t  =  0) 


(A p-  |p|)sign(p)|sin(7n/)| 
0 


if  (A p  —  |p|)sign(p)  sin(7n/)  >  0 
otherwise, 


(A p-  |p|)sign(p)|  sin(7rc/)| 
0 


if  (A p  —  |p|)sign(p)  sin(7n/)  <  0 
otherwise, 


where  these  are  orthonormalized  numerically  as  described  in  §5.2.  The  pressure  and  velocity 
for  all  modes  are  zero  everywhere,  initially. 

Stochastic  coefficient  initialization:  The  discrete  pdf  for  the  first  stochastic  coefficient  is 
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specified  with  four  samples  Yrl  =  [—0.18,  —  0.06,  0.04,  0.20]T,  (Fig.  [2]-c).  The  next  two 
coefficients  do  not  introduce  additional  uncertainty  because  we  specify  perfectly  correlated 
samples,  Yr> 2  =  Yr )3  =  e-  [—0.18,  —0.06,  0.04,  0.20]T,  where  e  is  a  small  constant  chosen  such 
that  Var(Y^j)  =  Var(yrii)  numerically.  This  means  that  the  inverse  of  the  covariance 
matrix  is  very  ill-conditioned  (or  numerically  singular),  so  the  pseudo- inverse  is  required 
during  the  initial  stages  of  the  simulation  (see  §5.1).  Also,  because  the  pdf  is  discrete,  we 
calculate  moments  using  the  biased  estimator  Cy^-.  ~  f  )T)r  YrtYrjl  instead  of  the  usual 
unbiased  estimator  C y.y.  «  ^2rYr  iYrj.  The  initial  helds  and  pdf  are  shown  in  Fig.  [2], 

6.1.2.  Lock-Exchange  Verification  Benchmark:  Results  and  Discussion 

The  outputs  from  the  stochastic  run  are  reported  in  Fig.  [3]  for  both  Grashoff  numbers. 
Comparisons  with  the  deterministic  runs  are  shown  in  Fig.  [4]  and  Fig.  [5]  for  Gr  =  4  x  104 
and  Gr  =  1.25  x  106,  respectively.  Finally,  the  evolution  of  the  differences  between  the 
stochastic  and  deterministic  runs  for  both  Grashoff  numbers  are  shown  in  Fig.  [6]. 

We  see  excellent  agreement  between  the  stochastic  realizations  and  the  deterministic 
runs  (Fig.  [5] -Fig.  [6])  for  this  challenging  benchmark.  Particularly,  for  the  lower  Grashoff 
number  flow,  the  local  error  is  less  than  0.2%  everywhere.  It  is  non-trivial  that  the  complex 
DO  implementation  is  capable  of  reproducing  multiple  deterministic  runs  in  a  single  simula¬ 
tion.  Thus,  based  on  these  results  and  many  other  tests  (not  shown),  we  conclude  that  our 
implementation  is  correct,  that  is,  we  are  solving  the  intended  equations. 

The  growth  of  differences  between  the  deterministic  and  stochastic  simulations  (Fig. 
[6])  should  be  explained  by  the  main  differences  between  the  stochastic  and  deterministic 
solvers,  in  particular  the  advection  schemes,  and  the  evolution  of  the  stochastic  coefficients. 
We  found  that  the  magnitude  of  the  differences  over  time  is  larger  for  coarser  space  and  time 
resolution  runs.  This  suggests  the  error  is  due  to  spatial  and/or  temporal  truncation  error. 
The  advection  scheme  does  not  contribute  significantly  to  the  error  at  the  reported  resolution, 
since  using  a  CDS  advection  scheme  instead  of  the  TVD*  scheme  for  the  stochastic  modes 
did  not  change  the  reported  results  significantly.  Reducing  the  time-step  to  At  =  1/1024 
reduced  the  error  at  the  final  time  from  ~2.1%  to  ~1.05%  for  the  higher  Grashoff  number 
flow.  This  indicates  that  the  error  is  dominated  by  a  temporal  truncation  compounded  error 
of  approximtely  O(At)  (as  expected,  see  §6.3).  The  primary  source  of  this  compounded  error 
may  be  from  the  evolution  of  the  stochastic  coefficients  and/or  the  modes. 

We  examine  the  errors  more  closely  to  determine  the  primary  source  of  the  small  numeri¬ 
cal  errors.  Note  that  the  density  differences  are  either  a  bit  too  small  (first  realization  in  Fig. 
[5])  or  a  bit  too  large  (last  three  realizations  in  Fig.  [5]),  causing  phase  errors  as  observed 
around  the  density  interface.  These  phase  errors  can  be  caused  by  small  errors  in  the  rel¬ 
ative  magnitudes  of  the  stochastic  coefficients.  We  found  that  different  orthonormalization 
strategies  (§5.2),  which  affect  the  relative  magnitudes  of  the  stochastic  coefficients,  can  also 
impact  the  magnitude  of  the  phase  errors.  In  particular,  when  using  a  Gram-Schmidt  or¬ 
thonormalization,  the  average  difference  was  greater  than  3%  (or  50%  worse)  for  the  high  Gr 
case  at  the  final  time.  Also,  considering  that  the  larger  Gr  number  flow  has  larger  stochastic 
coefficients  (Fig.  [3])  than  the  lower  Gr  flow,  the  same  relative  error  in  the  stochastic  coef¬ 
ficients  would  lead  to  a  larger  phase  error  (observed  when  comparing  Fig.  [4]  to  Fig.  [5]). 
Thus,  this  benchmark’s  results  indicate  that  the  primary  source  of  error  originates  from  the 
time  evolution  of  the  stochastic  coefficients. 
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Mode  2  Mode  3 


Streamlines 


Figure  3:  The  DO  mean,  modes  at  non-dimensional  time  t  =  5,  and  evolution  of  stochastic  coefficients, 
for  Gr  =  1.25  x  106  (top)  and  Gr  =  4  x  104  (bottom).  The  higher-Gr  flow  has  sharper  gradients,  and  its 
coefficients  are  larger  than  the  lower-Gr  flow.  Streamlines  shown  over  density  in  color.  Note  that  the  sign 
of  the  contribution  from  a  mode  depends  on  the  sign  of  the  stochastic  coefficients. 
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Final  Time:  Gr  =  1 .25e6  Final  Time:  Gr  =  4e4 


Figure  4:  The  four  DO  realizations  (top  row),  the  deterministic  runs  (middle  row),  and  the  DO  realizations 
minus  the  deterministic  runs  normalized  by  ||3>  deterministic^  (bottom  row)  for  Gr  =  4  x  104  (resolution 
256x256  with  512  x  5  time-steps).  The  stochastic  solver  result  agrees  with  the  deterministic  solver  results, 
which  validates  our  DO  numerical  schemes.  Streamlines  shown  over  density  in  color. 
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Figure  5:  As  Fig.  [4]  but  with  Gr  =  1.25  x  106.  The  higher-Gr  flow  with  sharper  gradients  has  larger  errors 
than  the  lower-Gr  flow,  but  they  are  still  small. 


Figure  6:  The  relative  errors  for  both  Gr  tend  to  grow  over  time  but  can  decrease.  In  both  cases,  the  DO 
mean  held  has  a  smaller  error  than  the  average  error  of  the  realizations. 
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In  summary,  using  the  proposed  benchmark,  we  have  verified  the  correct  implementation 
of  the  stochastic  code.  Also,  based  on  the  results,  we  suggest  that  the  accuracy  of  the  scheme 
will  benefit  most  from  improving  the  temporal  discretization  for  problems  that  require  long 
time  integration. 


6.2.  Effect  of  Advection  Scheme 

The  purpose  of  this  benchmark  is  to  test  the  three  different  advection  schemes  proposed 
(§4.1.2).  We  again  use  a  version  of  the  lock-exchange  problem  because  the  sharp  density 
interface  will  highlight  numerical  oscillations.  We  modify  the  problem  by  introducing  sym¬ 
metry,  which  should  be  maintained  numerically.  Finally,  for  simplicity  we  only  consider  a 
single  stochastic  mode  with  a  bimodal  continuous  pdf. 

General  setup:  The  Schmidt  number  is  kept  constant,  Sc  =  1,  and  we  present  results 
for  Gr  =  1.25  x  106.  Initially  the  velocity  is  zero  everywhere.  Free-slip  boundary  conditions 
are  used  at  the  boundaries  of  the  domain  (Fig.  [1]).  The  domain  is  discretized  using  Ax  = 
Ay  =  1/64,  and  At  =  1/256.  A  lower  resolution  is  used  here  compared  to  the  cases  in  §6.1 
in  order  to  highlight  the  symmetry  errors  and  numerical  oscillations 

Mean  initialization:  The  mean  density,  pressure,  and  velocity  are  all  zero  everywhere, 
initially. 

Mode  initialization:  The  density  profile  for  the  mode  is  the  same  normalized  hyperbolic 
tan  profile  used  in  §6.1.  The  pressure  and  velocity  for  this  mode  are  zero  everywhere,  initially. 

Stochastic  coefficient  initialization :  The  bimodal  Gaussian  continuous  pdf  is  represented 
by  10,000  samples.  To  ensure  that  any  asymmetry  comes  from  the  numerical  errors  only,  we 
exactly  enforce  symmetry  in  the  initial  conditions  as  follows.  We  first  generate  2,500  samples 
(Yr2 500)  from  a  zero  mean  Gaussian  distribution  with  standard  deviation  a  =  O.Ole1  ~  0.027. 
Next,  these  samples  are  duplicated  to  obtain  a  representation  of  the  bimodal  pdf  that  is 
exactly  symmetric  at  the  numerical  level: 


W,!  = 


Y _ —Y _ Y 

x  C2500  2  ’  r2500  2  ’  r2500 


-  —Y  +  - 

2  1  C2500  —  2 


These  are  the  10,000  samples  that  we  evolve.  They  are  illustrated  in  the  first  row  of  Fig. 
[7].  Of  course,  these  samples  are  correlated.  The  procedure  should  not  be  used  in  general; 
it  is  used  here  solely  to  evaluate  how  good  are  numerical  schemes  at  maintaining  symmetry. 


6.2.1.  Effect  of  Advection  Scheme:  Results 

The  result  of  this  simulation  is  shown  in  Fig.  [7]  for  the  three  advection  schemes.  While 
the  CDS  scheme  maintains  symmetry  of  the  mean,  modes,  pdf,  and  realizations,  some  clear 
numerical  oscillations  are  present,  particularly  evident  in  the  realizations.  While  there  are 
no  oscillations  in  the  TVD  scheme,  it  clearly  loses  symmetry  in  the  mean,  modes,  pdf,  and 
realizations,  as  can  be  seen  (with  aid  of  the  dashed  guide  lines)  in  Fig.  [7].  The  TVD* 
scheme  only  develops  minor  oscillations,  which  can  be  barely  detected  when  examining  the 
realizations,  and  completely  retains  symmetry.  Thus,  the  new  TVD*  scheme  is  the  preferred 
scheme  among  the  three. 

Initially,  we  can  represent  all  density  realizations  exactly  with  one  mode.  However,  as 
different  densities  evolve  at  different  rates,  one  mode  becomes  insufficient  to  represent  the 
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Figure  7:  Symmetric  lock-exchange  problem  (Gr  =  1.25  x  106)  with  grid  resolution  64  x  64  and  At  =  1/256 
using  various  advection  schemes  for  the  modes.  Only  the  stochastic  density  is  non-zero  initially,  with  a 
bi-modal  pdf(top  row).  The  two  most  extreme  realizations  (largest  and  smallest  stochastic  coefficients)  are 
plotted  in  each  case  (right  column).  Our  new  averaged  TVD*  scheme  only  has  minor  oscillations  and  retains 
symmetry  (third  row),  while  the  CDS  advection  scheme  suffers  from  large  oscillations  (white  dashed  circles, 
second  row),  and  the  one-sided  TVD  scheme  loses  symmetry  (red  dashed  circles,  last  row).  Streamlines 
shown  over  density  in  color. 
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uncertainty.  Hence,  spurious  gradients  can  appear  in  the  reconstructed  realizations.  We 
purposely  chose  to  use  only  one  mode  in  this  benchmark  to  also  illustrate  that  if  the  number 
of  modes  is  fixed,  errors  occur.  In  general,  we  do  not  keep  the  number  of  modes  fixed  [51]. 

In  summary,  we  found  that  the  TVD*  scheme  performs  adequately  and  that  our  DO 
implementation  can  reproduce  vastly  different  realizations  of  a  given  problem. 

6.3.  Numerical  Convergence  Analysis 

The  purpose  of  this  benchmark  is  mainly  to  show  that  the  implemented  scheme  is  converg¬ 
ing.  Here  we  use  the  classical  lid-driven  cavity  flow,  and  examine  the  numerical  convergence 
under  spatial  and  temporal  refinement  of  each  component  separately.  This  benchmark  does 
not  have  a  variable  density,  and  so  we  report  the  Reynolds  number  instead  of  the  Grashoff 
number.  We  also  completed  convergence  tests  with  density,  with  results  analogous  to  those 
presented  below. 

6.3.1.  Numerical  Convergence  Analysis:  Setup 


x-  (0,1) 


u  =  (1,0), 


x  =  (0,0)  x  -  (1,0) 


Figure  8:  The  lid-driven  cavity  flow  is  a  classical  benchmark  used  to  verify  convergence  of  numerical  imple¬ 
mentations.  Uncertainty  for  this  case  will  be  introduced  through  the  initial  conditions. 

General  setup:  We  present  results  for  a  Reynolds  number  of  Re  =  500,  although  other 
Reynolds  number  cases  (Re  G  [100, 1000])  were  also  studied,  giving  similar  results.  The  flow 
is  driven  by  a  deterministic  boundary  condition  at  the  top  of  the  enclosed  box  Fig.  [8],  with 
no-slip  velocity  boundary  conditions,  and  uncertain  initial  conditions.  The  finest  resolution 
uses  Ax  =  Ay  =  1/512,  and  At  =  1/4096,  which  is  sufficient  at  this  Reynolds  number 

[12,  11],  and  the  order  of  convergence  is  approximated  as  O  ~  log  /log(2). 

We  examine  the  difference  between  the  fine  and  coarse  space  resolutions  by  interpolating 
(using  splines)  the  fine  solution  onto  the  coarse  resolution  grid,  and  taking  the  L 2  norm  over 
the  interior  of  the  domain,  Vj  G  [0.25,0.75]  x  [0.25,0.75],  to  avoid  the  boundary  condition 
singularities  at  the  top  two  corners. 

Mean  initialization:  For  a  challenging  case,  the  mean  velocity  and  pressure  are  initially 
zero,  everywhere. 

Mode  initialization:  The  velocity  modes  are  initialized  by  specifying  the  stream  function 
iPm,n(x,  y )  =  Cm,n  sin(7nc)  sin(7rMa;)  sin(7n/)  sin(nNy), 
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where  Cm,n  =  yy  Ar^1  (f)^  ' '  +  iV^1  is  the  normalization  constant;  the  delta 

function  6(x)  takes  the  value  1  if  x  —  0  and  0  otherwise.  The  velocity  modes  are  then 
specified  as 


ILi 


Vi  ~  dx^{M,N)i  ’ 


where  (M,  N)  =  {(1, 1),  (1,2),  (1,  3)}.  The  initial  pressure  for  the  modes  is  specified  as  zero 
everywhere. 

Stochastic  coefficient  initialization :  The  pdf  is  created  using  5,000  samples  of  zero  mean 
Gaussian  distributions  with  variances  Var(Y)5000;j)  =  e^~Mi~Ni\  The  number  of  samples  is 
not  critical  in  this  case,  since  the  samples  are  the  same  from  one  run  to  the  next.  That  is,  we 
only  test  spatial  and  temporal  convergence,  not  stochastic  convergence.  Since  this  is,  again, 
a  numerical  test,  the  F*  samples  are  purposely  created  using  a  procedure  as  in  §6.2, 


r  Y  1 

y*  _  -‘*'5000,* 

Xr,i  ~  _y 

1  7*5000 

To  ensure  that  the  final  generated  samples  have  numerical  variances  exactly  as  specified,  we 
correct  the  samples  using  the  numerically  calculated  variance 

T  -  _  Y*  VV^Yr—) 

"  v'Vad'O;,)  ' 

where  Var9(arjj)  =  ^-j-  YH=i  ar,i  is  the  calculated  sample  variance.  The  initialization  for  this 
problem  can  be  seen  in  the  first  row  of  Fig.  [9]. 


6.3.2.  Numerical  Convergence  Analysis:  Results  and  Discussion 

Three  time  snapshots  of  the  reference  solution  are  shown  in  Fig.  [9].  While  the  marginal 
seems  to  remain  approximately  Gaussian,  the  sample-scatter  diagram  clearly  shows  the  very 
non-Gaussian  behavior  for  this  benchmark.  Also,  we  confirm  from  Fig.  [10]  that  the  variances 
of  the  modes  are  decreasing,  which  is  expected  since  this  problem  has  a  deterministic  steady- 
state  solution.  Thus,  the  stochastic  solution  is  behaving  as  expected. 

Next,  examining  the  convergence,  we  see  that  the  numerical  error  for  all  components  de¬ 
crease  with  temporal  and  spatial  refinement.  The  convergence  is  near  optimal  at  large  grid 
sizes  for  all  variables.  These  results  with  the  velocity  components  separated  are  tabulated 
in  Table  [1]  and  Table  [2],  which  also  shows  that  the  stochastic  coefficients  are  converging 
optimally.  Even  though  a  fourth-order  RK  method  is  used  to  advance  the  stochastic  coef¬ 
ficients,  the  total-DO  convergence  is  first  order  (based  on  choices  made  in  §3.5).  Thus,  we 
observe  near-optimal  convergence  for  all  variables,  which  suggests  that  the  implementation 
is  correct. 


6-4-  Stochastic  Convergence 

The  purpose  of  the  fourth  benchmark  is  to  further  assess  numerical  performance,  demon¬ 
strate  that  a  significant  number  of  modes  can  be  used,  and  study  the  effect  of  the  stochastic 
discretization.  Specifically,  we  quantify  the  effects  on  accuracy  of  the  number  of  stochastic 
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Figure  9:  Evolution  of  reference  solution  for  the  lid-driven  cavity  flow  (Re  =  500)  with  stochastic  initial 
conditions  over  time  (resolution  512x512  with  4096  x  5  time-steps).  This  case  demonstrates  that  the  DO 
representation  with  our  implementation  is  able  to  capture  and  evolve  non-Gaussian  statistics  (see  sample 
scatter  plot)  for  a  continuous  pdf.  The  marginal  pdfs  are  normalized  on  the  plot,  and  the  variance  can  be 
read  off  Fig.  [10].  Streamlines  shown  over  pressure  in  color.  The  sample  scatter  is  colored  by  the  s  =  3 
stochastic  coefficient. 


Time 


Figure  10:  Evolution  of  the  stochastic  energy  (Var(Y)))  for  the  lid-driven  cavity  reference  solution.  The 
tick-marks  on  the  time-axis  corresponds  to  the  time  snapshots  in  Fig.  [9]. 


Figure  11:  The  control  volume  size  is  held  fixed  at  Aa;  =  1/512  for  the  time  convergence  (left),  and  the  time 
step  size  is  held  fixed  at  At  =  1/4096  for  the  spatial  convergence  (right).  The  error  (||</>2iv  —  </>jv II2)  decreases 
with  both  temporal-  and  spatial-refinement  for  each  component,  and  convergence  is  near-optimal  (order  1 
in  time  and  2  in  space). 


Figure  12:  Laminar  vortex  shedding  over  a  square  cylinder  in  a  channel.  Here  uncertainty  originates  from  the 
initial  conditions  and  uncertain  vortex  shedding:  depending  on  the  perturbation,  the  first  vortex  could  either 
be  shed  above  or  below  the  cylinder.  Within  a  stochastic  framework,  however,  if  uncertainties  are  initially 
symmetric,  the  mean  and  modes  should  remain  symmetric,  and  this  is  evaluated  with  our  DO  numerics. 
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p 


u 


V 


Yi 


Nt 


Mean 


2048 

1024 

2048 

1024 

2048 

1024 

2048 

1024 


l|e||2 

2.9e-04 

5.7e-04 

2.7e-04 

5.4e-04 

2.5e-04 

5.0e-04 


O 

1.5 

1.0 

1.5 

1.0 

1.5 

1.0 


Mode  1 
||  ^  ||  2 
5.8e-04 
1.2e-03 
3.7e-04 
7.3e-04 
4.7e-04 
9.5e-04 
3.1e-04 
6.2e-04 


O 

1.5 

1.0 

1.5 

1.0 

1.5 

1.0 

1.5 

1.0 


Mode  2 
||  ^  ||  2 
3.0e-03 
5.9e-03 
1.6e-03 
3.1e-03 
1.8e-03 
3.5e-03 
5.1e-04 
1.0e-03 


O 

1.5 

0.99 

1.5 

.99 

1.5 

0.99 

1.5 

1.0 


Mode  3 
||  6 1|  2 
2.3e-03 
4.7e-03 
2.2e-03 
4.3e-03 
2.4e-03 
4.8e-03 
9.9e-04 
2.0e-03 


O 

1.5 

1.0 

1.5 

1.0 

1.5 

1.0 

1.5 

1.0 


Table  1:  Temporal  convergence  of  lid-driven  cavity  flow.  Tabulated  is  the  error  (e  =  \\4>2Nt  ~  4>Nt  H2)  between 
the  solutions  using  At  =  l/{2Nt)  and  using  At  =  l/iVt,  and  the  approximate  order  of  convergence  O.  The 
grid  size  is  fixed  at  Ax  =  1/512. 


samples  and  of  the  time-order  of  integration.  To  do  so,  we  extend  the  classic  shedding  of  vor¬ 
tices  by  a  uniform  flow  as  it  encounters  a  symmetric  obstacle  to  stochastic  DO  computations. 
Since  this  problem  is  symmetric,  the  stochastic  solution  can  only  lose  symmetry  if  numerical 
or  external  perturbations  initiate  the  non-symmetric  laminar  shedding  of  vortices.  However, 
if  perturbations  are  symmetric,  there  should  be  no  preferential  direction  for  vortex  shedding. 
Thus,  a  carefully  initialized  simulation  should  be  able  to  capture  symmetric  directions:  this 
provides  an  excellent  test  to  assess  numerical  performance. 

6-4- 1-  Stochastic  Convergence:  Setup 

The  benchmark  is  an  open  flow  in  a  frictionless  pipe  with  a  square  cylindrical  obstacle 
(Fig.  [12]),  which  is  a  classic  test  for  deterministic  flow  solvers. 

General  setup:  We  present  results  for  a  Reynolds  number  of  Re  =  100,  although  other 
cases  were  also  studied.  The  flow  is  driven  by  a  deterministic  uniform  inlet  boundary  con¬ 
dition  (left  of  domain),  with  slip  velocity  boundary  conditions  at  the  top  and  bottom,  open 
boundary  conditions  at  the  outlet,  and  symmetric  uncertain  initial  conditions.  All  simula¬ 
tions  use  a  resolution  of  336  x  63  in  space,  and  63  x  40  in  time.  We  choose  to  integrate 
until  t  —  40  because  this  allows  the  statistics  to  reach  steady  values.  At  t  =  40  the  mean 
velocity  has  traveled  through  the  domain  2.5  times. 

Mean  initialization:  The  mean  velocity  and  pressure  are  initially  zero,  everywhere. 

Mode  initialization:  The  exact  shape  of  the  initial  stochastic  perturbations  are  not  im¬ 
portant  since  they  are  advected  out  of  the  domain.  However,  to  maintain  symmetry,  per¬ 
turbations  have  to  be  symmetric.  We  initialize  the  velocity  modes  by  specifying  the  stream 
function 


'0m,7v(^,|/)  =  Cm,n  sin(nx/a)  sm^Mx/a)  sm(ny/b)  sm^Ny /b), 
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•  Nx 

Mean 
||  6 1|  2 

O 

Mode  1 
||e||2 

O 

Mode  2 

llel|2 

O 

Mode  3 
||  6 1|  2 

O 

P 

256 

1.8e-04 

2.2 

5.4e-04 

1.8 

1.9e-03 

2.3 

2.9e-03 

2.0 

128 

1.0e-03 

2.5 

2.1e-03 

2.0 

1.0e-02 

2.4 

l.le-02 

2.0 

64 

4.5e-03 

2.2 

9.6e-03 

2.2 

4.7e-02 

2.2 

3.0e-02 

1.5 

u 

256 

2.0e-04 

2.1 

5.2e-04 

2.1 

1.2e-03 

1.9 

3.6e-03 

2.1 

128 

1.0e-03 

2.4 

2.7e-03 

2.4 

4.9e-03 

2.1 

1.7e-02 

2.2 

64 

4.5e-03 

2.1 

1.3e-02 

2.3 

2.5e-02 

2.4 

6.7e-02 

2.0 

V 

256 

1.9e-04 

2.0 

5.6e-04 

2.1 

1.2e-03 

1.9 

4.0e-03 

2.1 

128 

8.6e-04 

2.2 

2.8e-03 

2.3 

5.1e-03 

2.0 

1.9e-02 

2.3 

64 

3.6e-03 

2.1 

1.4e-02 

2.4 

2.6e-02 

2.4 

7.8e-02 

2.0 

Yi 

256 

4.8e-04 

2.1 

7.2e-04 

2.1 

1.3e-03 

2.2 

128 

2.2e-03 

2.2 

3.4e-03 

2.2 

6.4e-03 

2.3 

64 

7.5e-03 

1.8 

1.5e-02 

2.1 

2.1e-02 

1.7 

Table  2:  Spatial  convergence  of  lid-driven  cavity  flow.  Tabulated  is  the  error  (e  =  ||^2JVX  —  </>jvJ|2)  between 
the  refined  (2 Nx  x  2NX)  and  present  (Nx  x  Nx )  grid,  and  the  approximate  order  of  convergence  O.  The  time 
step  is  fixed  at  its  smallest  value  At  =  1/4096. 


where  Cm,n  is  the  normalization  constant  (as  in  §6.3),  and  a  =  16,  b  =  3  are  the  width  and 
height  of  the  domain  respectively.  The  velocity  modes  are  then  specified  as 


d 

ui  = 


d 

Vi  = 


where 


(M,  N)  =  {(1, 1),  (2, 1),  (1,  2),  (3, 1),  (1,  3),  (2,  2),  (4, 1),  (1,4),  (3,  2),  (2,  3)}, 

and  Bm  is  a  smoothing  function  created  numerically  from  the  domain  mask.  BM  is  created 
from  the  mask  by  iteratively  averaging  each  control  volume  by  its  own  value  and  its  four 
neighbors  for  iterations.  That  is,  at  iteration  k 

=  t  (Bk„'(i,j)  +  £){,-'(!  -  1  ,j)  +  -  1) 

+Bm  l{i  +  1  ,j)  +  B^  1(i,j  +  1))  . 

The  initial  pressure  for  the  modes  is  specified  as  zero  everywhere. 

Stochastic  coefficient  initialization :  To  ensure  initial  symmetry,  the  samples  for  the 
stochastic  coefficients  are  created  using  the  same  procedure  as  in  §6.3,  using  the  variances 
Var(Yr>i )  =  e(2~Mi~Ni)  (note  the  difference  of  +1  in  the  exponential  from  §6.3).  The  reference 
solution  uses  105  samples  for  the  stochastic  coefficients,  and  a  4th  order  RK  time  integration 
scheme  (§3.5). 
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Figure  13:  Mean  field  and  first  5  modes  with  marginal  pdfs  at  final  non-dimensional  time  T=40  for  Re=100, 
and  the  evolution  of  the  mean,  (u,  u)^,  and  stochastic  energy,  Var(Fi),  for  the  reference  solution  (resolution 
63x336  with  63  x  40  time-steps).  Streamlines  shown  over  pressure  in  color.  Our  scheme  and  implementation 
retains  (anti)-symmetry  for  the  most  important  first  four  modes. 
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Energy  Energy 


Figure  14:  As  in  Fig.  [13],  but  showing  two  realizations  where  the  vortex  is  shed  in  opposite  directions.  The 
colorbar  for  the  pressure  is  as  on  Fig.  [13]. 


#  of  Samples 


Figure  15:  Energy  of  the  stochastic  coefficients,  Var(Y)),  over  time  for:  (top)  different  number  of  samples 
(0( 4)  time  integration);  and  (bottom),  different  time  integration  schemes  (10,000  samples).  Trends  are  well 
captured  in  all  cases,  but  there  are  noticeable  errors  for  less  energetic  modes  after  long  integration  times  when 
a  small  number  of  samples  and  a  low  order  time  discretization  scheme  is  used.  A  relatively  small  number 
of  samples  (10,000)  can  be  used  for  this  benchmark  since  nearly  identical  results  compared  to  100,000  were 
found  for  the  most  energetic  four  modes. 
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Samples:  1e2 
Order:  0(4) 


Mode  1 


Mode  2 


Mode  3 


Mode  4 


Mode  5 


1e3 

0(4) 


1e4 

0(4) 


1e4 

0(2) 


1e4 

0(1) 


1e5 

0(4) 


Figure  16:  Marginal  probability  density  functions  of  first  five  modes  for  the  flow  over  a  square  cylinder 
benchmark  at  final  time  (columns  1-3,  increasing  sample  sizes;  columns  4-5,  decreasing  order  in  time;  column 
6,  reference).  With  10,000  samples  for  the  stochastic  coefficients,  the  continuous  marginal  pdfs  are  well- 
represented,  although  the  bimodal  peaks  lose  some  symmetry.  The  marginals  have  similar  shapes  for  all 
sample  sizes,  but  the  best  representations  have  larger  sizes.  Overall,  the  order  of  the  time-integration 
scheme  has  less  effect  than  the  sample  size  (the  overall  order  in  time  is  one). 
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6-4-2.  Stochastic  Convergence:  Results  and  Discussion 

For  the  reference  solution,  we  find  that  excellent  symmetry  is  maintained  for  the  first  four 
stochastic  modes  (Fig.  [13]).  From  the  realizations  (Fig.  [14]),  the  scheme  clearly  captures 
both  shedding  directions.  We  find  that  for  fewer  samples  and  lower  time  integration  accuracy, 
the  symmetry  is  not  as  well  maintained  (not  shown).  However,  for  a  sufficient  number  of 
samples,  symmetry  is  maintained  for  all  Reynolds  numbers  we  tested.  From  Fig.  [13],  we 
see  that  variances  seem  to  reach  a  steady  value  after  an  initial  transient  period.  The  smaller 
variances  take  longer  to  reach  a  steady  value;  the  highest  modes  still  evolve  after  the  final 
time  step. 

For  our  sequential  implementation,  the  reference  DO  simulation  (Fig.  [13])  was  computed 
in  about  7.5  hours  on  a  2.4Ghz  computer.  Each  time-step  required  the  inversion  of  11 
pressure  equations,  as  opposed  to  the  111  required  for  the  scheme  without  the  stochastic 
pseudo  pressure.  This  roughly  translates  to  a  1000  %  increase  in  efficiency,  or  about  3  days 
saving  in  terms  of  computational  time  for  this  computer.  Since  then,  we  have  used  as  many 
as  s  =  25  modes  for  data  assimilation  applications  [56,  57]. 

Sample  sizes:  Examining  the  evolution  of  the  stochastic  coefficients  for  different  sample 
sizes  (top  of  Fig.  [15]),  we  see  that  the  differences  are  larger  for  higher  modes.  For  the 
most  energetic  first  four  coefficients,  the  results  agree  well  with  the  reference  solution  when 
using  >  1,000  samples.  After  the  5th  coefficient,  the  differences  become  larger,  and  in 
the  10th  coefficient,  the  differences  begin  sooner  and  have  larger  relative  magnitudes.  Note 
that  the  logarithmic  scale  magnifies  smaller  errors,  and  the  differences  for  the  10th  are  fully 
insignificant  when  compared  to  the  variance  of  the  first  mode.  Thus,  with  a  relatively  small 
number  of  samples,  it  is  possible  to  capture  the  stochastic  coefficient’s  variances  accurately. 
Also,  it  appears  as  though  the  solution  converges  when  the  number  of  samples  increases, 
which  increases  our  confidence  in  this  approach. 

Not  only  are  the  evolution  of  the  variances  well-represented  by  a  smaller  number  of 
samples,  the  shape  of  the  marginals  are  also  well-reproduced  (Fig.  [16]).  We  observe  that  the 
marginals  are  well-reproduced  using  >  1,  000  samples,  and  for  all  time  integration  schemes 
(Fig.  [16]).  The  magnitudes  of  the  bimodal  peaks  in  the  first  mode  are  more  symmetric 
with  increased  sample  sizes,  which  suggests  using  upwards  of  10,000  samples  is,  perhaps, 
advised.  In  fact,  using  such  a  number  of  samples  does  not  affect  the  overall  cost  of  the 
solution  scheme.  That  the  marginals  are  well-represented  using  smaller  sizes  is  encouraging, 
since  it  suggests  that  a  small  number  of  samples  (that  won’t  negatively  impact  the  efficiency 
of  the  scheme)  could  be  used  for  problems  with  large  s:  this  would  not  substantially  affect 
accuracy. 

Order  in  time:  We  are  here  only  varying  the  order  of  the  ODE  solver  used  to  evolve  the 
stochastic  coefficients:  the  overall  order  of  the  DO  scheme  is  kept  fixed  at  first  order  (see 
§3.4).  Examining  results  (bottom  of  Fig.  [15]),  we  see  that  the  second  and  fourth  order 
ODE  solver  agree  well  for  all  coefficients.  A  second  order  solver  thus  seems  sufficient. 

In  summary,  we  simulated  a  stochastic  version  of  the  classical  flow  over  a  square  cylinder 
benchmark.  We  found  that  our  DO  solver  maintains  excellent  symmetry,  and  the  use  of 
pseudo-pressures  reduces  the  cost  by  1000%  when  using  10  stochastic  modes.  We  also  found 
that  the  statistics  is  converging  with  the  number  of  samples  and  that  accuracy  benefits 
from  using  a  second-order  solver  for  the  coefficients,  even  though  the  overall  DO  scheme  is 
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limited  to  first  order.  Finally,  smaller  sample  sizes  for  the  coefficients  still  produced  accurate 
time-evolving  statistic  and  marginal  pdfs,  which  is  encouraging. 

7.  Conclusions  and  Future  Research 

We  have  derived  efficient  computational  schemes  for  the  DO  methodology  applied  to 
unsteady  stochastic  Navier-Stokes  and  Boussinesq  equations,  and  illustrated  and  studied 
the  numerical  aspects  of  these  schemes.  For  the  discretizations  in  time,  we  developed  semi- 
implicit  projection  methods  for  the  mean  and  the  modes,  and  we  employed  time-marching 
schemes  of  first  to  fourth  order  for  the  stochastic  coefficients.  For  the  discretizations  of  the 
physical  space,  we  employ  conservative  second-order  finite- volumes,  with  a  special  treatment 
for  the  advection  terms  based  on  a  TVD  scheme  with  monotonized  central  symmetric  flux 
limiter.  Several  alternate  discretizations  in  time  and  space  were  outlined  and  illustrated  by 
comparisons.  We  also  addressed  several  numerical  issues  specific  to  the  DO  method  for  fluid 
and  ocean  flows.  In  particular  we  have  shown:  how  to  define  pseudo-stochastic  pressures 
to  reduce  the  number  of  matrix  inversions  for  the  pressure  from  0(s2)  to  O(s);  how  to 
treat  advection  by  the  stochastic  modes  using  symmetric  approximate  TVD  schemes;  how 
to  deal  with  singular  subspace  covariances  by  generalized  inversion;  and  how  to  maintain 
orthonormal  modes  at  the  numerical  level  and  so  account  for  truncation  and  round-off  errors 
during  time  integration.  Finally,  we  evaluated  our  schemes  using  a  varied  set  of  stochastic 
flows,  hence  illustrating  robustness  but  also  providing  benchmarks  for  future  schemes  and 
implementations. 

Using  an  asymmetric  Dirac-stochastic  lock-exchange  benchmark,  we  found  excellent  agree¬ 
ment  among  multiple  realizations  from  a  deterministic  code  and  the  realizations  generated 
from  a  single  stochastic  simulation.  This  validated  our  numerical  implementation,  and  con¬ 
firmed  that  the  pseudo-stochastic  pressure  approach  works  in  practice.  Using  a  symmetric 
version  of  the  lock-exchange  benchmark,  we  showed  that  the  symmetric  TVD  advection 
scheme  (TVD*)  works  well,  while  the  CDS  scheme  suffers  from  numerical  oscillations,  and 
the  non-symmetric  TVD  scheme  loses  symmetry. 

Using  a  lid-driven  cavity  flow  with  uncertain  initial  conditions,  we  showed  that  each  com¬ 
ponent  converged  near-optimally  under  both  time  and  space  refinement.  Additionally,  we 
showed  that  even  when  a  higher-order  time  integrator  is  used  for  the  stochastic  coefficients, 
their  accuracy  is  still  limited  to  approximately  first  order,  as  discussed  in  §3.4.  This  bench¬ 
mark  also  demonstrated  an  expected  decay  in  the  variance  of  the  stochastic  coefficients, 
which  were  shown  to  be  very  non-Gaussian. 

Finally,  using  a  stochastic  flow  past  a  square  cylinder  in  a  confined  channel,  we  showed 
that  the  stochastic  coefficients  converged  with  increased  samples  sizes  and  orders  of  time  in¬ 
tegration.  We  found  that  using  >  10,  000  samples  and  a  second  order  time-integrator  yielded 
adequate  performance.  This  benchmark  also  demonstrated  that  our  discretized  DO  method¬ 
ology  successfully  captures  both  vortex  shedding  directions,  resulting  in  a  fully  symmetric 
mean  held.  Also,  we  noted  that  using  our  newly  defined  pseudo-stochastic  pressure,  we  were 
1000%  more  efficient  computationally  when  using  10  modes. 

Possible  future  studies  include  the  efficient  extension  of  our  present  framework  to  hows 
with  uncertain  parameters  in  the  governing  equations  and  with  stochastic  forcing  at  the 
boundary  and  in  the  interior.  Also,  our  results  suggest  that  simulations  with  long  time 
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integration  will  benefit  from  more  accurate  time-integration  scheme.  Then,  while  our  pro¬ 
posed  symmetric  TVD-based  advection  and  orthonormalization  schemes  for  the  stochastic 
modes  were  shown  to  perform  adequately,  more  optimal  treatments  are  possible.  Alternative 
means  of  discretizing  the  stochastic  coefficients  are  worth  exploring  for  applications  where 
rare  events  are  important.  Finally,  specific  multi-resolution  DO  schemes  can  be  derived  for 
multiscale  stochastic  fluid  and  ocean  flows.  Unstructured  grids  as  well  as  nested  approaches 
[8,  18]  will  then  be  useful.  The  utilization  of  the  present  schemes  as  well  as  these  future  ad¬ 
vances  will  allow  the  prediction  of  uncertainty  in  realistic  simulations  of  multi-physics  fluid 
systems,  over  a  wide  range  of  applications,  from  micro-nano  fluid  engineering  to  multiscale 
ocean  and  climate  studies. 
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Appendix  A.  DO  Navier-Stokes  Equations 

In  this  Appendix  we  present  the  stochastic  DO  equations  for  Boussinesq  dynamics,  and 
briefly  summarize  the  required  steps  (from  [50,  49,  52])  to  obtain  these  equations. 

Starting  with  a  generalized  Karhunen-Loeve  expansion  [32,  49],  one  decomposes  the  so¬ 
lution  of  eqns.  (l)-(3)  into  a  Dynamically  Orthogonal  (DO)  field  expansion  [50,  49,  51]  for 
the  velocity,  density  and  pressure2  fields 


$(x,t;o;)  =  $(x,£)  + 

1=1 

&  =  $  +  Yi&i, 

s  s  s 

P(x,  t;  u)  =  p(x,  t)  +  y,  Yi(t,  c u)pi(x,  t)  +  EE  Yi(t;  u)Yj(t\ uj)pij(x,  t), 

i= 1  2=1  j= 1 

p  =  p  +  YjPi  +  YjYjPij, 

where  $  =  [cf*1,  $2, . .  ,]7  =  [u,  p]T  is  the  vector  of  prognostic  state  variables.  The  scalar 
s  =  s(t)  defines  the  time-dependent  dimension  of  the  stochastic  subspace,  i.e.  s  is  a  discrete 
number  of  stochastic  terms  retained  from  a  complete  expansion.  The  field  functions  <&j(x,  t) 
are  the  s  orthonormal  deterministic  modes  and  the  cj)  are  their  s  zero- mean  stochastic 


(A.l) 

(A.2) 


2For  convenience,  we  changed  the  sign  of  from  the  original  definition  in  the  Sapsis  and  Lermusiaux 
[50] 
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coefficients,  in  general  non-Gaussian.  In  our  notation,  we  use  the  Einstein  summation  ex¬ 
clusively  for  summations  related  to  the  stochastic  expansion.  The  decomposition  of  pressure 
into  a  mean,  linear  modal  and  quadratic  modal  component  follows  from  the  Pressure  Poisson 
Equations  (see  §3.2  or  Sapsis  and  Lermusiaux  [50]).  Since  both  the  modes  and  stochastic 
coefficients  are  functions  of  time,  a  redundancy  arises,  which  is  resolved  by  the  DO  condition, 


<I>; 


d*i 

dt 


=  0,  Vi,  j  G  [1,  2, ... ,  s], 


v 


(A.3) 


where  the  inner  product  is  defined  as  (a,  b)^  =  4D  for  arbitrary  vectors  of 

spatial  functions  a  =  [a1,  a2, . .  .]T  and  b  =  [61,  b2, . .  ,]T. 

Using  the  DO  condition,  an  exact  set  of  equations  can  be  obtained  that  governs  the 
evolution  of  the  mean,  modes,  and  stochastic  coefficients  of  the  generalized  Karhunen-Loeve 
expansion.  The  only  approximation  arises  from  the  truncation  of  the  DO  expansion  to  s(t) 
terms.  First  substitute  (A.l)  into  eqns.  (l)-(3)  to  obtain  (using  a  Langevin  notation): 


<9u 

~dt 

dp 

dt 


+ 

+ 


dY>  _u  v 

-7Tui  +  Y* 

at 

dYi 

~MPi  +  Y‘ 


dui 

dt 

dpi 

dt 


Cu  (u  +  Y^,  p  +  Y^,  p,  x,  t]  u)  , 
Cp  ( p  +  Y^,  u  +  FjUj,  x,  t]  u)  , 


(A.4) 


and  DO  decomposed  versions  of  (2)-(3).  It  is  from  these  equations,  within  which  the  DO 
decomposition  was  inserted,  that  the  equations  for  the  mean,  modes  and  their  coefficients  are 
obtained,  using  the  expectation  operator,  the  spatial  inner  product,  and  eqns.  (A.2)-(A.3). 

Mean.  To  obtain  a  rate  of  change  for  the  mean  fields,  the  idea  is  to  eliminate  the  random 
components  in  the  left-hand-sides  of  eqns.  (A.4).  Hence,  taking  the  expectation  of  eqns.  (A.4) 
it  can  be  found  that  the  evolution  of  the  mean  fields  are  governed  by 


V  ■  u  =  0,  xfD, 

—  -  — = V2u  =  -V  ■  (uu)  -  Vp  +  pe9 
dt  v  Gr 

-  C Y-Yj  (Vpg  +  V  •  (UjU;))  ,  x  G  V, 

~  ScVCdV2f)  =  -V  ’  ^  _  Cv‘Y)V  '  (UjPi^  x  G 
boundary  conditions  given  by 


dp 

dt 

with  the  initial  and 


u  (x,  0) 
P  (x,  0) 
u 
<9u 
dn 
P 
dp 
dn 


Uo, 

xeP, 

do, 

xGD, 

g  Di 

x  G  dVD, 

Sn, 

x  G  dVN, 

9  Dp, 

x  G  dVDp 

9np, 

x  G  dVNp 

(A.5) 


(A. 6) 


(A.7) 
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where  C^y.  =  Eu  [YjYj]  is  an  element  of  the  covariance  matrix  in  the  stochastic/error 
subspace.  Deterministic  initial  and  boundary  conditions  ( •  quantities)  are  assigned  to 
the  mean.  Note  that  in  general  the  vector  V  ■  (u.jUj)  differs  from  V  ■  (UjU,-)  (recall  that 


{UjUi)q  = 


9u, 


7^ 


duirUjq '  ~  _  dvjUi  /  dvi 


,  OVjUi  / 

e.g.,  -±-  7^ 


dx  r  '  dxr  ’  dy  '  dy  ' 

Modes.  The  evolution  of  the  modes  is  also  obtained  from  eqns.  (A. 4).  To  do  so,  the 
idea  is  to  eliminate  the  random  coefficients  in  front  of  the  time  derivatives  of  the  modes. 
The  essential  steps  are  to  multiply  these  equations  with  a  stochastic  coefficient  Yj,  apply 
the  expectation  operator,  and  substitute  an  expression  for  Eu  [^-Y)] ,  which  is  obtained  by 
projecting  the  equation  unto  and  imposing  the  DO  condition.  From  this  the  following 
governing  evolution  equations  for  the  modes  can  be  found: 


V  ■  u,:  =  0,  x  e  Z>, 

(9u ' 

Q_j_  =  Q i  ~  (Qi,  *&j) v  ui  ’  x  ^  (A. 8) 

^  =  Q?  -  (Qi,  ®j)v  py  x  e  v, 


where 


Qi  =  [Q",Q f]  =  [c y^^Y^  ,  cv'./a  /yoy_ 

Q"  =  -7=  V2u i  -  V  ■  (ujii)  -  V  ■  (uuj)  -  Vpi  +  p^9 


VG~r 

'Y,Y,MYiYmYn  (Vp„m  +  V  •  (Unll,m)  )  , 


cr1 


QP  =  _ _ 

^  Sc\fGr 


V2pi  -  V  •  (Ujp)  -  V  ■  (u Pi 


~  (V  '  ( UnPm ) ) , 


and  My.ymyn  =  Eu  [Y]YmYn]  is  a  third  moment.  The  right-hand-sides  in  eqns.  (A. 8)  corre¬ 
spond  the  total  rate  of  change  of  the  subspace  (without  a  DO  condition)  minus  the  projection 
of  this  rate  of  change  on  the  subspace  itself  (which  is  subtracted  to  ensure  the  DO  condition). 
We  note  that  in  general  V  •  (uUj)  7^  V  ■  (Ujli)  since  7^  for  example. 

The  initial  and  boundary  conditions  for  the  modes  are  obtained  from  those  of  the  full 
stochastic  fields,  eqns.  (2)  and  (3),  but  reduced  to  their  dominant  initial  error  subspace  of 
size  so, 


Uj  (x,  0)  —  Ujj0  (x) ,  xeD, 
pi  (x,  0)  =  Pi, 0  (x) ,  x  G  V, 


(A.9) 
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and 


U  i  =  gi,D, 


dui 

dn 


gi,JV, 


Pi  “  .9i,Dp, 
dpt  _ 
dn  ~  9i'Np' 


x  G  &Dd, 

X  G  <9Xhv, 
x  G  <9X>up, 
x  G  <9£>jV 


(A.10) 


Coefficients.  Finally,  to  obtain  the  evolution  of  the  stochastic  coefficients  from  eqns.  (A. 4), 
the  idea  is  to  eliminate  the  modes  in  the  term  containing  the  time  derivatives  of  the  ran¬ 
dom  coefficients.  To  do  so,  project  the  evolution  eqns.  (A. 4)  onto  each  mode  i,  apply  the 
DO  conditions,  and  essentially  impose  that  each  coefficient  is  of  zero  mean.  The  resulting 
governing  stochastic  ODEs  are 


<m 

dt 


Ym  (V pnm  T  V  •  (unUm),  U f) ^  (YmYn  C}'myn) 
—  (V  ■  (u nPm),  Pi)rp  (frnYn  ~  Cymyn)  , 


(A.ll) 


where  £  =  [£u,  £fT  and  Fm  =  [F"  ,  Ff]T  with 

1 


Fu  = 

m  VOr 

i 

Fp  = 


V2u  m 


V  ■  (umu)  —  V  ■  uu„ 


Vpm  +  Pm?9, 


ScVG r 


V2Pm  ~  V  •  (Ump)  -  V  ■  (Upm), 


The  initial  conditions  for  the  coefficients  are  obtained  from  those  of  the  full  stochastic 
holds,  eqn.  (2),  by  projection  onto  each  initial  mode  i,  and  removal  of  the  mean.  This 
leads  to: 


Yi  {to]  u)  -  ($0  -  $0  ,  ,  ,A  12x 

=  (uo  —  u0  ,  U^o)^  +  (po  —  Po  ,  Pi, o)  ,  WGfl. 

In  this  Appendix  we  stated  the  DO  equations  for  the  mean,  modes,  and  stochastic  coef¬ 
ficients  of  the  stochastic  Boussinesq  equations. 

Appendix  B.  Numerical  orthonomalization  procedure 

In  this  Appendix,  we  present  our  newly  developed,  numerically  efficient,  orthonormaliza¬ 
tion  procedure  based  on  the  discussion  in  §5.2. 

Consider  the  pictorial  analogy  showing  different  strategies  when  orthonormalizing  the  ba¬ 
sis  (Fig.  [B. IT] ) .  In  case  (a),  the  argument  is  that  the  coefficients  were  evolved  correctly,  and 
all  the  error  lies  in  the  basis  evolution.  This  is  not  true  because  we  evolve  the  coefficients  by 
taking  an  inner  product  of  the  equations  using  the  bases  at  the  old  time-step.  In  case  (b),  the 
argument  is  that  the  energy  and  ‘direction’  of  the  realizations  were  correctly  evolved.  When 
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a) 

Orthonormalize 
basis  without 
modifying 
coefficients 


b) 

Orthonormalize 
basis  and  retain 
realization 
energy  and 
'direction' 


|Cov(V)=Cov(X)| 
(Y,Y) 

. . 

®  (Y'.Y) 


Figure  B.17:  There  are  different  ways  to  treat  the  stochastic  coefficients  when  orthononnalizing  the  modes. 
In  this  pictorial  analogy,  the  standard  2D  cartesian  bases  and  coordinates  represent  the  inodes  and  stochastic 
coefficients,  respectively. 


correcting  the  bases,  the  direction  and  energy  are  maintained,  but  the  actual  realization  is 
modified.  In  case  (c),  the  argument  is  that  the  realizations  were  evolved  correctly.  However, 
it  is  unlikely  that  the  errors  in  the  modes  and  coefficients  cancel  during  the  time  evolution. 
Also,  while  not  evident  from  the  pictorial  analogy,  the  energy  of  the  coefficients  change  in 
case  (c).  This  is  because  the  length  of  the  basis  changes  when  orthonormalized,  and  this 
length  change  is  then  accounted  for  in  the  coefficients  in  order  to  maintain  the  coordinate 
of  the  realization.  The  case  (b)  strategy  is  a  compromise  between  cases  (a)  and  (c),  and  it 
results  in  the  lowest  error  for  the  lock-exchange  cases  in  §6.1.  The  errors  were  the  largest 
for  case  (a),  and  the  errors  for  case  (c)  were  marginally  larger  compared  to  case  (b). 

To  be  clear,  case  (b)  is  shown  to  perform  best  for  the  specific  cases  in  §6.1,  and  these 
results  are  not  proven  to  be  universally  applicable.  Also,  we  made  no  attempt  to  provide  an 
exhaustive  list  of  all  the  orthonormalization  strategies.  For  example,  arguing  that  the  coef¬ 
ficients  still  live  in  the  bases  of  the  old  timestep,  one  can  consider  projecting  the  coefficients 
from  this  old  basis  onto  the  new  basis  before  orthonormalization.  Therefore,  the  issue  of  the 
optimal  orthonormalization  strategy  is  left  open  to  future  research.  Nonetheless,  we  believe 
our  current  approach  is  sufficient  for  most  applications  since  no  strategy  will  influence  the 
scheme’s  overall  order  of  accuracy.  Any  improved  strategy  will  only  reduce  the  size  of  the 
error,  which  will  be  of  order  At  smaller  than  the  truncation  error.  Therefore,  any  accuracy 
issues  due  to  orthonormalization  can  be  resolved  by  reducing  the  size  of  the  timestep  when 
using  the  current  scheme.  Our  orthonormalization  procedure  for  case  (b)  is  described  next. 

For  case  (b),  we  want  the  new  orthonormalized  basis  to:  approximately  recover  the  spe¬ 
cific  realizations  before  the  orthonormalization;  and  exactly  maintain  the  stochastic  energy. 
Mathematically  we,  ideally,  want  both: 


rO\T 


Tr(Cy*  Y  )  =  Tr(C‘*’° 

\  lrl±r  n/  \  lrlIr 


(B.l) 

(B.2) 


where  the  orthonormalized  basis  is  defined  such  that  —  & iji  and  the  trace 

operator  Tr  simply  sums  the  diagonal  entries  of  a  matrix.  In  order  to  preserve  the  second 
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property,  we  first  rotate  the  stochastic  system  such  that  the  covariance  matrix  is  diagonal: 


C  kYiYj  =  E  [YriiYrJ]  =  VdcDdcV£c,  (B.3) 

Yrf  =  Y„Vnc,  (B.4) 

K?  =  *h,iVDC,  (B.5) 

which  gives  the  decorrclated  (or  DC)  stochastic  coefficients.  By  starting  with  a  decorrelated 
system,  we  have  a  unique  and  convenient  reference  frame  which  we  can  use  to  enforce  the 
Tr(Cf?r)  property. 

Next  we  perform  the  orthonomalization  as  follows 


=  <#£f ,  Sgf  >B  =  VmDmVtm 

y°c  =  y°cvmV oy, 

d^. 


(B.6) 

(B.7) 

(B.8) 


The  above  could  be  accomplished  by  directly  calculating  the  Singular  Value  Decomposition 
of  the  matrix,  but  the  above  procedure  has  the  advantage  of  enabling  the  user  to  specify 
any  desired  (possibly  non-linear)  inner  product.  In  either  case,  the  two  approaches  are  similar 
in  efficiency  and  give  the  same  results.  To  verify  that  (B.l)  is  satisfied,  we  substitute  the 
above  expression  in  (B.l) 


qPc(ydc)t 

h,i  V  r,z  ) 


(y°cVm-/d„)T  , 
«°fv«v/iD’\/D^vii(r®c)T, 

qDC(ydc)t 

h,i  v  r,i  )  i 


since  the  eigen- vectors  are  orthonormal  =  I.  Therefore  the  first  property,  (B.l),  is 

exactly  satisfied  at  this  stage.  Stopping  here  would  give  the  case  (c)  strategy. 

Finally,  we  decorrelate  the  samples  again,  while  ensuring  that  the  stochastic  energy  is 
preserved: 


<d;?f  =  Eiuyuy]  =  VoDoVj, 


rOCvOC 


Y°  =  Yoc\, 


CM 


K  =  vo- 


'Tr(D 


DC 


TV(Do) 


(B.9) 

(B.10) 

(B.ll) 


Note  that  Tr(Do)  =  Ti'(Dz)cDm),  but  that  Do  ^  DocDm-  This  ensures  that  (B.2)  is 
maintained,  but  now  we  can  only  approximately  satisfy  (B.l).  The  over-all  correction  is 
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then  as  follows: 


Y°  =  RiVDCVMV/D^Voy 

ITt(Ddc) 

(B.  12) 

Tr(D0)  ’ 

Ki  =  ^mV^Vma/d-Vo- 

(B.13) 

This  procedure  is  robust  even  if  c  is  singular,  however  it  does  require  a  non-singular 
Mg  which  will  normally  be  the  case  if  modes  are  properly  initialized.  Also,  the  number  of 
floating  point  operations  is  dominated  by  the  calculation  of  Mg  if  the  number  of  spatial  points 
exceeds  the  number  of  realizations.  Nonetheless,  the  overall  cost  of  the  orthonormalization 
is  inexpensive  compared  to  the  inversion  of  the  pressure  Poisson  equations. 
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